{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampling the potential energy surface "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction \n",
    "\n",
    "This interactive notebook demonstrates how to utilize the Potential Energy Surface (PES) samplers algorithm of qiskit chemistry to generate the dissociation profile of a molecule. We use the Born-Oppenhemier Potential Energy Surface (BOPES)and demonstrate how to exploit bootstrapping and extrapolation to reduce the total number of function evaluations in computing the PES using the Variational Quantum Eigensolver (VQE). \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import common packages\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import partial\n",
    "\n",
    "# qiskit\n",
    "from qiskit.aqua import QuantumInstance\n",
    "from qiskit import BasicAer\n",
    "from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE, IQPE\n",
    "from qiskit.aqua.components.optimizers import SLSQP\n",
    "from qiskit.circuit.library import ExcitationPreserving\n",
    "from qiskit import BasicAer\n",
    "from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE\n",
    "from qiskit.aqua.components.optimizers import SLSQP\n",
    "\n",
    "# chemistry components\n",
    "from qiskit.chemistry.components.initial_states import HartreeFock\n",
    "from qiskit.chemistry.components.variational_forms import UCCSD\n",
    "from qiskit.chemistry.drivers import PySCFDriver, UnitsType\n",
    "from qiskit.chemistry.core import Hamiltonian, TransformationType, QubitMappingType\n",
    "from qiskit.aqua.algorithms import VQAlgorithm, VQE, MinimumEigensolver\n",
    "from qiskit.chemistry.transformations import FermionicTransformation\n",
    "from qiskit.chemistry.drivers import PySCFDriver\n",
    "from qiskit.chemistry.algorithms.ground_state_solvers import GroundStateEigensolver\n",
    "from qiskit.chemistry.algorithms.pes_samplers.bopes_sampler import BOPESSampler\n",
    "from qiskit.chemistry.drivers.molecule import Molecule\n",
    "from qiskit.chemistry.algorithms.pes_samplers.extrapolator import *\n",
    "\n",
    "import warnings\n",
    "warnings.simplefilter('ignore', np.RankWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we use the H2 molecule as a model system for testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ft = FermionicTransformation()\n",
    "driver = PySCFDriver()\n",
    "solver = VQE(quantum_instance=\n",
    "             QuantumInstance(backend=BasicAer.get_backend('statevector_simulator')))\n",
    "me_gsc = GroundStateEigensolver(ft, solver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "stretch1 = partial(Molecule.absolute_stretching, atom_pair=(1, 0))\n",
    "mol = Molecule(geometry=[('H', [0., 0., 0.]),\n",
    "                        ('H', [0., 0., 0.3])],\n",
    "                       degrees_of_freedom=[stretch1],\n",
    "                       )\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = PySCFDriver(molecule=mol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.3])]\n"
     ]
    }
   ],
   "source": [
    "print(mol.geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a perturbation to the molecule along the absolute_stretching dof"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.5])]\n",
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.8999999999999999])]\n"
     ]
    }
   ],
   "source": [
    "mol.perturbations = [0.2]\n",
    "print(mol.geometry)\n",
    "\n",
    "mol.perturbations = [0.6]\n",
    "print(mol.geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate bond dissociation profile using BOPES Sampler\n",
    "\n",
    "Here, we pass the molecular information and the VQE to a built-in type called the BOPES Sampler. The BOPES Sampler allows the computation of the potential energy surface for a specified set of degrees of freedom/points of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First we compare no bootstrapping vs bootstrapping \n",
    "\n",
    "Bootstrapping the BOPES sampler involves utilizing the optimal variational parameters for a given degree of freedom, say r (ex. interatomic distance) as the initial point for VQE at a later degree of freedom, say r + $\\epsilon$. By default, if bootstrapping is set to True, all previous optimal parameters are used as initial points for the next runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance1 = partial(Molecule.absolute_distance, atom_pair=(1, 0))\n",
    "mol = Molecule(geometry=[('H', [0., 0., 0.]),\n",
    "                        ('H', [0., 0., 0.3])],\n",
    "                       degrees_of_freedom=[distance1],\n",
    "                       )\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = PySCFDriver(molecule=mol)\n",
    "\n",
    "# Specify degree of freedom (points of interest)\n",
    "points = np.linspace(0.25, 2, 30)\n",
    "results_full = {} # full dictionary of results for each condition\n",
    "results = {} # dictionary of (point,energy) results for each condition\n",
    "conditions = {False: 'no bootstrapping', True: 'bootstrapping'}\n",
    "\n",
    "\n",
    "for value, bootstrap in conditions.items():\n",
    "    # define instance to sampler\n",
    "    bs = BOPESSampler(\n",
    "        gss=me_gsc\n",
    "        ,bootstrap=value\n",
    "        ,num_bootstrap=None\n",
    "        ,extrapolator=None)\n",
    "    # execute\n",
    "    res = bs.sample(driver,points)\n",
    "    results_full[f'{bootstrap}'] =  res.raw_results\n",
    "    results[f'points_{bootstrap}'] = res.points\n",
    "    results[f'energies_{bootstrap}'] = res.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare to classical eigensolver\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define numpy solver\n",
    "solver_numpy = NumPyMinimumEigensolver()\n",
    "me_gsc_numpy = GroundStateEigensolver(ft, solver_numpy)\n",
    "bs_classical = BOPESSampler(\n",
    "               gss=me_gsc_numpy\n",
    "               ,bootstrap=False\n",
    "               ,num_bootstrap=None\n",
    "               ,extrapolator=None)\n",
    "# execute\n",
    "res_np = bs_classical.sample(driver, points)\n",
    "results_full['np'] =  res_np.raw_results\n",
    "results['points_np'] = res_np.points\n",
    "results['energies_np'] = res_np.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5nklEQVR4nO3deXxU9b3/8ddnZrKwhJCNJWwJkBCWQICACKIgaLUVEbfqlQpyXdC2Fvuz1V7789a2v9atat3riraKO25FLSCICIphDzuEACEsWQgEQkhm5vv7Y07iECZhCJmcLJ/n4zEPzpxz5pz3nAzzmXO+53yPGGNQSimlauOwO4BSSqmmTQuFUkqpOmmhUEopVSctFEoppeqkhUIppVSdtFAopZSqkxYK1ayIyPMi8n9tWO9nIjKtnq/tKSJHRcTZ0Lkamoj8WUQKRWR/zdwislhEbrY7o2p8otdRqKZCRHKBzoAb8AAbgdeBF4wxXhujnRHrfdxsjFlgd5YzISI9gS1AL2PMwQDTFwP/Msa81NjZlL10j0I1NZOMMVFAL+BB4B7gZXsjtQwi4jrNLD2BokBFQrVuWihUk2SMOWyM+Rj4KTBNRAYBiMhsEfmzNRwvIp+KSImIFIvI1yLisKbdIyJ7RaRURLaIyARrfISIPCEi+dbjCRGJqFqviEwWkTUickREdojIJdb46sMuItJHRL4UkSLrMM0bItLRmvZPfF+4n1iHbX4rIkkiYqq+qEUkUUQ+tjJvF5Fb/Nb/BxF5R0Ret7JvEJHM2raTtdw7RSTHyvKI3zaYLiLfiMjjIlIE/EFEoq1lF4jILhH5vYg4RGQiMB9ItHLPrpk7wLpniMgmETkkIl+ISK/6/K1V06eFQjVpxpgVQB4wNsDk/2NNS8B3yOp/ACMi/YBfACOsvZMfAbnWa+4DRgEZwBBgJPB7ABEZie9Q12+AjsD5fq/zJ8BfgUSgP9AD+IOV92fAbnx7Ru2NMQ8HeP1bVu5E4GrgLyJyod/0y615OgIfA08HWIa/KUAmMAyYDMzwm3YOkINv+/w/4CkgGugNXADcCNxkHSa7FMi3ck+va4UiMhnf9r4S3/b/GphzmpyqmdJCoZqDfCA2wPhKoCu+Y+qVxpivja/RzQNEAANEJMwYk2uM2WG95gbgj8aYg8aYAuAB4GfWtP8GXjHGzDfGeI0xe40xm2uu1Biz3ZrnhLWMx/B96Z6WiPQAxgD3GGPKjTFrgJfwfWFXWWqMmWeM8QD/xFfQ6vKQMabYGLMbeAK43m9avjHmKWOMG6gArgN+Z4wpNcbkAn/ze/9nYibwV2PMJmvZfwEydK+iZdJCoZqDbkBxgPGPANuB/1iHXu4F3xc5MAvfr/yDIvKWiCRar0kEdvktY5c1Dnx7Bjs4DRHpbC1zr4gcAf4FxAf5XhKBYmNMaY0M3fye7/cbLgMiT9O+sKfGshJrmRYPhHHq+/dfd7B6AX+3DvuV4Pv7SD2XpZo4LRSqSROREfi+fJbWnGb9Kv4/xpje+A7X/LqqLcIY86Yx5jx8X2gGeMh6Wb41rkpPaxz4vlT7BBHrL9Yy040xHYCp+L4kq6PV8dp8IFZEompk2BvEemvTo8ay8v2e+2cpxLcXVvP912fde4DbjDEd/R5tjDHL6rEs1cRpoVBNkoh0EJHL8B2r/5cxZn2AeS4Tkb4iIsBhfIecvCLST0QutBqpy4HjQNXptXOA34tIgojEA/fj2yMA39lVN4nIBKuBt5uIpAWIFwUcBQ6LSDd8bRr+DuBrAziFMWYPsAz4q4hEishgfIe8/hVo/iD9RkRirMNavwLermXdHuAd4P+JSJR1mOjX9Vz388DvRGQggNVIfk394qumTguFamo+EZFSfL9Y78N3/P+mWuZNARbg+9JeDjxrjFmEr33iQXy/oPcDnYDfWa/5M5AFrAPWA6uscVUN5zcBj+MrPF9x8q/vKg/gazg+DPwb+KDG9L/iK0YlInJ3gNdfDyTh++U/F/jfs7zm4iNgJbDGylPX6cS/BI7ha+BeCrwJvHKmKzTGzMW3l/aWdfgtG19juGqB9II7pZoxETFAitUuo1RI6B6FUkqpOmmhUEopVSc99KSUUqpOukehlFKqTqfrJKzZiY+PN0lJSXbHUEqpZmXlypWFxpiEQNNaXKFISkoiKyvL7hhKKdWsiMiu2qbpoSellFJ10kKhlFKqTloolFJK1anFtVEopc5OZWUleXl5lJeX2x1FhUBkZCTdu3cnLCws6NdooVBKnSQvL4+oqCiSkpLw9beoWgpjDEVFReTl5ZGcnBz06/TQk1LqJOXl5cTFxWmRaIFEhLi4uDPeW9RCoZQ6hRaJlqs+f1stFJa8g7nc9+oU5i193e4oSinVpGihsDhE+NixnRU75tkdRSl1FqZPn85777131stZs2YN8+bV/X2Qm5vLm2++edbrCsbzzz/P66/b80NWC4UlMaEXHTxeCisO2B1FKdUEnG2hcLvdDZpn5syZ3HjjjQ26zGBpofCT4HZyyHvE7hhKtWq5ubn079+fW265hYEDB3LxxRdz/PhxwPflPWrUKAYPHsyUKVM4dOhQwGUsWLCAzMxMUlNT+fTTTwFfI/1NN91Eeno6Q4cOZdGiRbWOr6io4P777+ftt98mIyODt99+m6+++oqMjAwyMjIYOnQopaWl3HvvvXz99ddkZGTw+OOPM3v2bC6//HIuvPBCJkyYwNGjR5kwYQLDhg0jPT2djz76qPo9pqWlccMNN9C/f3+uvvpqysrKAF83RL/97W9JT09n5MiRbN/uuyfVH/7wBx599FEAxo0bxz333MPIkSNJTU3l66+/BqCsrIxrr72WAQMGMGXKFM4555wG6dJIT4/1E2vastd51O4YSjUZD3yygY35DfvjaUBiB/530sA659m2bRtz5szhxRdf5Nprr+X9999n6tSp3HjjjTz11FNccMEF3H///TzwwAM88cQTp7w+NzeXFStWsGPHDsaPH8/27dt55plnEBHWr1/P5s2bufjii9m6dWut4//4xz+SlZXF008/DcCkSZN45plnGDNmDEePHiUyMpIHH3yQRx99tLoYzZ49m1WrVrFu3TpiY2Nxu93MnTuXDh06UFhYyKhRo7j88ssB2LJlCy+//DJjxoxhxowZPPvss9x9t+/OudHR0axfv57XX3+dWbNmVS/fn9vtZsWKFcybN48HHniABQsW8OyzzxITE8PGjRvJzs4mIyPjLP5SP9A9Cj+xrhgKXFBRccLuKEq1asnJydVfcsOHDyc3N5fDhw9TUlLCBRdcAMC0adNYsmRJwNdfe+21OBwOUlJS6N27N5s3b2bp0qVMnToVgLS0NHr16sXWrVtrHV/TmDFj+PWvf82TTz5JSUkJLlfg39kXXXQRsbGxgO+6hf/5n/9h8ODBTJw4kb1793LggO/wdo8ePRgzZgwAU6dOZenSpdXLuP7666v/Xb58ecD1XHnllSdtH4ClS5dy3XXXATBo0CAGDx4c8LVnSvco/HRq053Kyjw271rN4JRRdsdRynan++UfKhEREdXDTqez+tBTsGqeAtoQp/vee++9/OQnP2HevHmMGTOGL774IuB87dq1qx5+4403KCgoYOXKlYSFhZGUlFR9DUNdGWsb9le1jZxOZ4O3h9SkexR+usWmALB190qbkyilaoqOjiYmJqb6ePw///nP6r2Lmt599128Xi87duwgJyeHfv36MXbsWN544w0Atm7dyu7du+scHxUVRWlpafUyd+zYQXp6Ovfccw8jRoxg8+bNp8xT0+HDh+nUqRNhYWEsWrSIXbt+6Ml79+7d1XsLb775Juedd171tLfffrv633PPPTfobTRmzBjeeecdADZu3Mj69euDfm1ddI/CT0q3YXDgn+wu3GR3FKVUAK+99hozZ86krKyM3r178+qrrwacr2fPnowcOZIjR47w/PPPExkZyR133MHtt99Oeno6LpeL2bNnExERUev48ePH8+CDD5KRkcHvfvc7li5dyqJFi3A4HAwcOJBLL70Uh8OB0+lkyJAhTJ8+nZiYmJNy3HDDDUyaNIn09HQyMzNJS0urntavXz+eeeYZZsyYwYABA7j99turpx06dIjBgwcTERHBnDlzgt4+d9xxB9OmTWPAgAGkpaUxcOBAoqOjz3Arn6rF3TM7MzPT1LeV/1hZKee+cy4/MX34600fNXAypZqHTZs20b9/f7tjtGi5ublcdtllZGdnnzKt6uZr8fHxZ7xcj8dDZWUlkZGR7Nixg4kTJ7JlyxbCw8NPmi/Q31hEVhpjMgMtV/co/LRrG0WCx1DkKbQ7ilJKnbGysjLGjx9PZWUlxhieffbZU4pEfWihqCHOE0YxeoqsUip0kpKSAu5NANVnMNVHVFRUSG4FrY3ZNcQSRaEztGcQKKVUc6KFooa4sASKXA5KSvXwk1JKgRaKU3Rp3wuA7O3f2pxEKaWaBi0UNfRM8J0JsCN/rc1JlFKqadBCUUP/5JEA5JdstzmJUq1Xbm4ugwYNOuvlPPHEE9Wd7dVm9uzZ5Ofnn/W6gjF69OhGWU9D00JRQ9/u6UR4DQXHG+eDo5QKnbMtFB6Pp0HzLFu2rEGX11i0UNTgcDrp5IZiT+Dui5VSjcPtdp/SDffChQsZOnQo6enpzJgxgxMnfB14Bhr/5JNPkp+fz/jx4xk/fjwej4fp06czaNAg0tPTefzxx3nvvffIysrihhtuICMjg+PHj5OUlMQ999zDsGHDePfdd3nxxRcZMWIEQ4YM4aqrrqouPNOnT2fmzJmndGc+e/ZsJk+ezLhx40hJSeGBBx6ofk/t27cHYPHixYwbN46rr766urvxqouf582bR1paGsOHD+fOO+/ksssua8zNHpAt11GISCzwNpAE5ALXGmMCfjOLSAdgI/ChMeYXjZEvzhtJkePMOiFTqkX67F7Y3zD9BVXrkg6XPnja2Wp2w/3YY4/xj3/8g4ULF5KamsqNN97Ic889x8yZM5k+ffop42fNmsVjjz3GokWLiI+PZ+XKlezdu7f6+oWSkhI6duzI008/zaOPPkpm5g8XJcfFxbFq1SoAioqKuOWWWwD4/e9/z8svv8wvf/lLIHB35gArVqwgOzubtm3bMmLECH7yk5+ctHyA1atXs2HDBhITExkzZgzffPMNmZmZ3HbbbSxZsoTk5OTqXmTtZtcexb3AQmNMCrDQel6bPwGB+xIOkVhHNAUuL94G3u1USgWvZjfcCxcuJDk5mdTUVOCHbsa3bNkScHxNvXv3Jicnh1/+8pd8/vnndOjQodZ1//SnP60ezs7OZuzYsaSnp/PGG2+wYcOG6mmBujMHX1fjcXFxtGnThiuvvPKkLsSrjBw5ku7du+NwOMjIyCA3N5fNmzfTu3dvkpOTAZpMobDryuzJwDhr+DVgMXBPzZlEZDjQGfgcCNgHSSjER3blmPcgew7soFdiamOtVqmmJ4hf/qFSs3vtjh07UlRUVO/lxcTEsHbtWr744guef/553nnnHV555ZWA8/p3FT59+nQ+/PBDhgwZwuzZs1m8eHGtGaueB9PNec2u1EPdVfjZsGuPorMxZp81vB9fMTiJiDiAvwF3n25hInKriGSJSFZBQcFZh+sa3RuAjTu/O+tlKaXqp2Y33JmZmeTm5lYf3qnqZrxfv34BxwMndQNeWFiI1+vlqquu4s9//nP1oaXTdRVeWlpK165dqaysrO6OvEqg7swB5s+fT3FxMcePH+fDDz+s3jM6nX79+pGTk1PdjUdVd+N2C9kehYgsALoEmHSf/xNjjBGRQF3Y3gHMM8bkne6mI8aYF4AXwNd7bP0S/6B3l8FwaC65BwL3xaKUCr2a3XA/+eSTjBo1imuuuQa3282IESOYOXMmERERvPrqq6eMB7j11lu55JJLSExM5IknnuCmm27C6/UC8Ne//hX4oVG6TZs2Ae8m96c//YlzzjmHhIQEzjnnnJOKSqDuzMF3WOmqq64iLy+PqVOnntI+UZs2bdrw7LPPcskll9CuXTtGjBhxVtuwodjSzbiIbAHGGWP2iUhXYLExpl+Ned4AxgJeoD0QDjxrjKmrPeOsuhmvcqBoLxM/vYQrGcAD05pGRVeqsWg348GZPn06l112GVdfffVJ42fPnn3SvbbP1NGjR2nfvj3GGH7+85+TkpLCXXfd1RCRq51pN+N2HXr6GJhmDU8DTrn5gzHmBmNMT2NMEr7DT6+frkg0lM5x3Yj2eCk8caAxVqeUUtVefPFFMjIyGDhwIIcPH+a2226zO5JtexRxwDtAT2AXvtNji0UkE5hpjLm5xvzTgcxgTo9tiD0KgCkvDKaNcfHmbavOellKNSe6R9HyNYsbFxljioAJAcZnATcHGD8bmB3yYH5iacceZ+0NXEop1Vroldm1iHXGUuCCiooTdkdRSilbaaGoRae23XGLsHHn93ZHUUopW2mhqEW3WN+Fdlt3axuFUqp100JRi9TuwwHIK95icxKllLKXFopaDOgzAocxHDy2x+4oSillKy0UtWgb2Y4Et6HIrffOVqqx5ebm0r9/f2655RYGDhzIxRdfzPHjxxk3bhxVp78XFhaSlJQE+C5yu+KKK7joootISkri6aef5rHHHmPo0KGMGjWK4uJiAMaNG8evfvUrMjIyGDRoECtWrMDr9ZKSkkJV9z9er5e+ffvSEN0BtRR2dQrYLMR7wynmmN0xlLLNQyseYnPx5gZdZlpsGveMPKUP0FNs27aNOXPm8OKLL3Lttdfy/vvv1zl/dnY2q1evpry8nL59+/LQQw+xevVq7rrrLl5//XVmzZoFQFlZGWvWrGHJkiXMmDGD7Oxspk6dyhtvvMGsWbNYsGABQ4YMISEhoSHebougexR1iCWKQmfT7dFRqZYsOTmZjIwMAIYPH17dUV5txo8fT1RUFAkJCURHRzNp0iQA0tPTT3ptVdfd559/PkeOHKGkpIQZM2bw+uuvA/DKK69w0003Nfj7ac50j6IOseEJFMshDh0uICZaf12o1ieYX/6hUrMb7uPHj+Nyuao79SsvL691fofDUf3c4XCc1IV3oC7Ae/ToQefOnfnyyy9ZsWLFKb3Etna6R1GHLlG+m4dk7zi1R0mlVONLSkpi5cqVALz33nv1WkZV191Lly4lOjqa6OhoAG6++WamTp3KNddcg9PpbJjALYQWijr0TPD1hbI9f63NSZRSAHfffTfPPfccQ4cOpbCwfieaREZGMnToUGbOnMnLL79cPf7yyy/n6NGjetgpAFs6BQylhuoUEGD77mymLLqe65zDuW/q7AZZplJNXUvuFHDcuHGn3B+7SlZWFnfddRdff/21DckaV7PoFLC56N2tP5FeQ2Flvt1RlFIh9OCDD/Lcc89p20QttFDUweF00sktFHsP2R1FKdUA/O937e/ee+/l3nsb5XY3zZK2UZxGnImkyFF++hmVakFa2iFp9YP6/G21UJxGrERT4PLi9XjsjqJUo4iMjKSoqEiLRQtkjKGoqKj63t7B0kNPpxEfmUiZ9wC79m8juVua3XGUCrnu3buTl5enXVi0UJGRkXTv3v2MXqOF4jQSO/aB4tVs3PmtFgrVKoSFhZGcnGx3DNWE6KGn0+jddTAAuw5utDmJUkrZQwvFaQzscw4A+0t32pxEKaXsoYeeTiMhJpEYj5ciz0G7oyillC20UAQh3u2kmCN2x1BKKVvooacgxJp2FDkq7Y6hlFK20EIRhNiwOA66oPxEmd1RlFKq0WmhCELntj3wiLBx50q7oyilVKPTQhGEbrGpAGzd3TC90iqlVHOihSIIqT19Pe/uLd5qcxKllGp8WiiCMCB5OE5jOFC2x+4oSinV6PT02CBERrSlkxuKPUV2R1FKqUanexRBivOGUSzH7I6hlFKNTgtFkGKlA4Uu7WpcKdX6aKEIUlx4Jw45HRSV7Lc7ilJKNSotFEHqEuXrdnn99uU2J1FKqcalhSJIvToNACBn33qbkyilVOPSQhGk/r1GAJBfst3mJEop1bj09NggJSWm0dbrpbByn91RlFKqUWmhCJLD6STB7aDYlNgdRSmlGpUeejoDcd5IiqXc7hhKKdWotFCcgVhnDAdcBq9Hr6dQSrUethQKEYkVkfkiss36N6aW+XqKyH9EZJOIbBSRpEaOepL4yETKHULO3k12xlBKqUZl1x7FvcBCY0wKsNB6HsjrwCPGmP7ASMDWG1cnduwDwMbc7+yMoZRSjcquQjEZeM0afg24ouYMIjIAcBlj5gMYY44aY2y9xVzfxCEA7C7QPQqlVOthV6HobIypOs90P9A5wDypQImIfCAiq0XkERFxBlqYiNwqIlkiklVQUBCqzAzqc64vcOnOkK1DKaWampCdHisiC4AuASbd5//EGGNExASYzwWMBYYCu4G3genAyzVnNMa8ALwAkJmZGWhZDSImOoFYt5diT+iKkVJKNTUhKxTGmIm1TRORAyLS1RizT0S6ErjtIQ9YY4zJsV7zITCKAIWiMcV7XBRTamcEpZRqVHYdevoYmGYNTwM+CjDP90BHEUmwnl8IbGyEbHWKpR2Fjgq7YyilVKOxq1A8CFwkItuAidZzRCRTRF4CMMZ4gLuBhSKyHhDgRZvyVotzxVPoEsrK9SZGSqnWwZYuPIwxRcCEAOOzgJv9ns8HBjditNPq1K4HnvKdbNzxPZkDx9kdRymlQk6vzD5D3WP7AbB97yqbkyilVOPQQnGG0nplArCnaIvNSZRSqnFooThDaUnDcRlDQVme3VGUUqpRaDfjZyg8PIJObijyFtsdRSmlGoXuUdRDnCecQ+hZT0qp1kELRT3EOqLZ7/Jod+NKqVZBC0U9JEelUep0sGLjQrujKKVUyGmhqIfMlEsA+G7zv21OopRSoaeFoh7GDPkx7bxetpesszuKUkqFXFCFQkT+JiIDQx2muXC5wkiujGA3hXZHUUqpkAt2j2IT8IKIfCciM0UkOpShmoNeru7sCjMUley3O4pSSoVUUIXCGPOSMWYMcCOQBKwTkTdFZHwowzVlaZ1H4BFh0cr37Y6ilFIhFXQbhXV3uTTrUQisBX4tIm+FKFuTdkHGNQBk531tcxKllAqtoK7MFpHHgUnAQuAvxpgV1qSHRKRVdnqU3C2NxEpDrifH7ihKKRVSwXbhsQ74vTEm0OXIIxswT7PSyxvNFlcJXo8HhzPg7byVUqrZC/bQ01qgn4gM83v0ERGXMeZwKAM2Zb3b9aPY5WDd9uV2R1FKqZAJtlA8C3wLvIDvLnPLgXeBLSJycYiyNXlDe/tuC75sw8c2J1FKqdAJtlDkA0ONMZnGmOHAUCAHuAh4OFThmrqxQycT4TVsK15jdxSllAqZYAtFqjFmQ9UTY8xGIM0Y06pbcttGtiO50sVu7wG7oyilVMgEWyg2ishzInKB9XjWGhcBVIYwX5PXy9GVneEeSo+V2B1FKaVCIthCMQ3YDsyyHjnAdHxFotVedAeQmpBJpV54p5RqwU5bKKwL7eYZY/5mjJliPR41xpQZY7zGmKONkLPJGpt+BQDrdi22NYdSSoXKaQuFMcYDeLV/p8D69x5OJ7eXnce32x1FKaVCItgL7o4C60VkPvxwD1BjzJ0hSdXMJLmj2OU8YncMpZQKiWALxQfWQwWQ1LYvK7xr2bJzNf2Sh9odRymlGlRQhcIY85qItAF6GmNaZd9OdRnccxzv5K7l6/VztVAopVqcYG9cNAlYA3xuPc8QEb0c2TJu+JW4jGFzwfd2R1FKqQYX7Omxf8DX+V8JgDFmDdA7JImaoej2sSRXONjt2Wd3FKWUanDBForKAJ3/eRs6THPW09GZnDA35SfK7I6ilFINKthCsUFE/gtwikiKiDwFLAthrmYnJTaDEw5hyaqP7I6ilFINKthC8UtgIHACmAMcwXeFtrKcO+ByAFbtXGBzEqWUaljBnvVUBtxnPVQAGSmjifnGy84TelKYUqplCfZWqKnA3UCS/2uMMReGJlbz43A6SXa3JddZYncUpZRqUMFecPcu8DzwEuAJXZzmrVdkMqvYxK78rfRKTLU7jlJKNYhg2yjcxpjnjDErjDErqx4hTdYMpXc7H4DFq9+1OYlSSjWcYAvFJyJyh4h0FZHYqkdIkzVDF2ZejcMYNh/4zu4oSinVYII99DTN+vc3fuMMetHdSeI6dqFXpbDLm2d3FKWUajDBnvWUHOogLUVP4skKP4jbXYnLFWZ3HKWUOmt1HnoSkd/6DV9TY9pfQhWqOesbPYhjDgfL131mdxSllGoQp2ujuM5v+Hc1pl1S35VabRzzRWSb9W9MLfM9LCIbRGSTiDwpIlLfdTaWkf1+DMCKrVoolFItw+kKhdQyHOj5mbgXWGiMSQEWWs9PXrjIaGAMMBgYBIwALjiLdTaKkQMnEuXxsuPIRrujKKVUgzhdoTC1DAd6fiYmA69Zw68BV9Sy7kggHIgAwoADZ7HORuFyhZFcGcluKbY7ilJKNYjTFYohInJEREqBwdZw1fP0s1hvZ2NMVZ/c+4HONWcwxiwHFgH7rMcXxphNgRYmIreKSJaIZBUUFJxFrIaRFN6T3WGG/YV77I6ilFJnrc5CYYxxGmM6GGOijDEua7jqeZ2n9IjIAhHJDvCYXGMdhgB7JyLSF+gPdAe6AReKyNhacr5gjMk0xmQmJCSc5i2H3oCu52JEWLxKL7xTSjV/wV5wd8aMMRONMYMCPD4CDohIVwDr34MBFjEF+NYYc9QYcxT4DDg3VHkb0gXDrgUgO/8bm5MopdTZC1mhOI2P+eEivmlAoJs47AYuEBGXiITha8gOeOipqeneKYkeFYZdJ3bZHUUppc6aXYXiQeAiEdkGTLSeIyKZIvKSNc97wA5gPbAWWGuM+cSOsPXR08SwM+w4Xo/2oaiUat6C7cKjQRljioAJAcZnATdbwx7gtkaO1mB6R/Xnm4rlZG1cxMj0iXbHUUqperNrj6LFy+xzMQDfbv7U5iRKKXV2tFCEyHkZk2jr9bK9ZJ3dUZRS6qzYcuipNQgPjyC5MpzdFNodRSmlzoruUYRQT1c3doV7OXTY/osAlVKqvrRQhFBawkjcIny5Ui+8U0o1X1ooQmj8sGsRY8jK/dzuKEopVW9aKEIouVsaAyvCWOnN0esplFLNlhaKEBsRPZp9YcK/v3nt9DMrpVQTpIUixK678LeEew0LtrxpdxSllKoXLRQhlpjQiyGV7Vjt3Ef5iTK74yil1BnTQtEIRnf5EYecDt5d+KTdUZRS6oxpoWgEP53wa6I8Xr7OazZ9GiqlVDUtFI0gql1HhnriWRNWQlHJfrvjKKXUGdFC0UjGJV/DcYeDOQsftTuKUkqdES0UjWTKuNtIcHv5rnCx3VGUUuqMaKFoJC5XGMPoSXZEOTv3brY7jlJKBU0LRSP6cfp/4xbh7cUP2x1FKaWCpoWiEV048mp6VUDWsZV2R1FKqaBpoWhkw8L7syXCy6qNX9kdRSmlgqKFopFdOepOAD5YoRffKaWaBy0UjSyj33mknXCyqnKL9iirlGoWtFDYYET7EewJFxaseMfuKEopdVpaKGzw03G/wWUMX2yYbXcUpZQ6LS0UNuiVmEp6RRtWSh4VFSfsjqOUUnXSQmGTUfHjKXI5+GDxs3ZHUUqpOmmhsMn1E35DW6+Xr3I/sDuKUkrVSQuFTWKiE8iojGF1WDGHjxbbHUcppWqlhcJG5/e8nGMOB28t0B5llVJNlxYKG101/hfEuL0sP7DA7ihKqWauouIEBYfyQ7JsV0iWqoISGdGWYSaRJRH72LM/hx5detsdSSnVBFRUnGD3/m3kHdzKgeJcDh07yOHjhZSeOMQxTynHvMcoMycokwqOObyUOgxHHUJqhZP3bl3b4Hm0UNjs4rSfsXDHI7y9+GHuvu55u+MopULE7a5k974t7Ni7gfyi7RSW5nH4RCGllSWUeo5yTMopFTdHnF4OOwSvyKkLcUBbvEQB7YyDdt4wYr0RtDVtaCdRdI1JDkl2LRQ2u+TcG/j75of5/sR3dkdRStVT6bESNu9cSc6+9RwoyaWobB+HK4o47C3liJRT4nRT7BTcAb782zq9RCNEeZ3Ee9vQy7Slg6sDHSJi6di2K3Htu9Ippgdd45JI7JRMVLuOjf7+tFDYzOF0MtyZwqeubWRv/45Bfc+xO5JSqob8gl1s3vk9uw9sZN+RnRQd388hzyFKpIwip4dip2BqFIF2Li8xHqGjN4w+no4MI5qOEZ2Ib9+NLjHJ9OiUSnLiAGKi4216V8HTQtEEXJH5cz5ZdRfvL3ucQX3fsjuOUq2O1+Nhe142G3OWk1uwgQNHd1FQWUCRHOOgy8MR58nn/YS5DAlArDeM/p5oYhzxxLfpRmJMH3p16U/fHoNJiEm0582EgBaKJmBk+kRSlgsrZANudyUuV5jdkZRqkfILdpG1aQE79q0m/2gORe4iCh1l7Hd5Oe74oRg4HIZOTkjwhDPUE0eCqzOd2veie3w/UnsMpU+PQa3q/6kWiiZiXMyFvFi2kGfm3s2vrvm73XGUara8Hg9bdq1h7fav2Fmwjn1luznoPcR+VwVFrh+Kgctp6GIg3hNJb3cMCW260SM2ldQemaT3OZd2baNsfBdNixaKJuKOKY+w6JXhzHUvYGrJfuI6drE7klJN3v7CPXy7fh6b960g71gO+0wxe8PcHPPbO2jr8pJY6aSfJ5YurkR6xg5gQK9RDO03lsiItjambz60UDQRLlcY01N/we93PcWjc2/jrzd9ZHckpZoMt7uSrE2LWLtjMbmHNpBfkU++8zj7w35oQG7v9NLDHcYId2e6tu1F74TBpPc5n/5JGTicThvTN39aKJqQyeNu5ZMXXmV+2A6u37qMwamj7Y6kVKNzuytZsWEBK7f+h5zDG8jzHmRXWGV1G4JDDIlO6OFtzyjTg97xgxnebyKDeo/UghAithQKEbkG+APQHxhpjMmqZb5LgL8DTuAlY8yDjRbSJneO+xszvrmVJxf9H15KXW53HKVCyuvxsHrzElZs+YwdxevJ8+4nN6yi+tBRhNOQ5HVyrrsrvTqk0r/7uZw7+FI6RjX9U0pbErv2KLKBK4F/1DaDiDiBZ4CLgDzgexH52BizsXEi2mNw6mguWtqbTyN38tHiF5g87la7IynVYPYX7uHLlW+xIX8ZuRW72Bl2glLr1NMwlyGpwsFId2eSOvRnaO8LOXfwpdqO0ATYUiiMMZsAJNAl6j8YCWw3xuRY874FTAZadKEA+PWU51n+wUW8tvVpfnLeTa3qNDzVcng9HrI2fsmyTZ+w/fA6dlHErjCDEUHE0NMhDPXEk9QmjSFJ4xgz5DI906iJasptFN2APX7P84CAly2LyK3ArQA9e/YMfbIQS4hJ5IqoCbx8fJGeLquajbLyY8z/bg6rcuezszyHnWFllFh7C+1cXvpURHKZ6cmgrqO5cPh1dInvYXNiFayQFQoRWQAEOsfzPmNMg57SY4x5AXgBIDMz0zTksu3yiyv/xld6uqxqwkqPlbBgxVuszJ3PjsqdbA+roNzhO0rQzWkY6I6ld5v+jEy5lPMyLtM942YsZIXCGDPxLBexF/D/ydHdGtcq6Omyqqk5fLSYL5b/izV5i8ipzGV7WCUnHAIO6CUw2t2ZAfEjmTD8Bvr2HGR3XNWAmvKhp++BFBFJxlcgrgP+y95IjUtPl1V2qqg4wfwVc1i+/RO2VeSwPaySCoevfSFJhLGeRAbGjeKizKn0Sky1O64KIbtOj50CPAUkAP8WkTXGmB+JSCK+02B/bIxxi8gvgC/wnR77ijFmgx157aSny6rG4vV4WLFxIYvXv83m0vVsDTvmOyNJIFngAk83BsWP4aIRU/UmW62MXWc9zQXmBhifD/zY7/k8YF4jRmty9HRZFUo5ezbw2fezyS78lq3OYg5afSF1cnnJcMcxqOMoLhkxjd49BtqcVNlJjGkRbb/VMjMzTVZWwOv3mq2CQ/lc88FFxHqcvDNjpTYKqnorP1HGF8v/xfKcT9nszmVHhO//f5THS7/K9vSLGsSF6deTOWC8XuXcyojISmNMZqBpTbmNQln0dFl1NrbuWsOn377EhpIsNoWVUup04HAYUnAy2ZvCuSmTuGjk9YSHR9gdVTVRWiiaCT1dVgWr/EQZ//5mNt/t/DebvbvZGe4bH+fyMtgdy5COo7ls9G3azqCCpoWimdDTZVVd9uzbxoffPMfa4uVsdB2h1OnA6TSkelxMIYWx/a5mwoir9XCSqhctFM2I/+myP92ylIx+59kdSdnE6/GwbN1nLFz/LzaUb2ZruBuPCDHWXsPQ2LFMGnMbiQm97I6qWgAtFM3MneP+xi3f3MJ9S27n2aj39fz1VuRYWSkfff0CK/Z8zkb2sc+6F0OSwI+8yZzXdwqXjv6ZnuygGpwWimZmcOpofr3rFh7Oe4lZn17N81d9Rue4bnbHUiGyZ38OHy59mjVFy8gOL6XM4SDcaehfEcmFbTP48YibGZwyyu6YqoXTQtEM/fSiWRz75DBPFr3Lne9fxj+un6/987cgqzZ/zbzvXyT7+Hq2hFfitg4pDa+MY3jXCVxx/kziojvbHVO1IloomqkZk/6XY++W8CLzuXPOpTx/42LaRrazO5aqB6/HwxffvslXW95lg2cnudZZSj3EMNHTi7EpV/Lj0TfqISVlG73grpl78I0ZvOH+nlEnonhuxlf6ZdJMHCsrZe5Xz/Jd3udkOw5S6HLgMIZ+FS4GRQ7g4qE3MSr9IrtjqlZEL7hrwe694RXKXruGuRGbuevVi/n7jAV6CmQTlV+wiw+WPMXqom/YEHaEYw4HkS7DwIp2TI46h8ljfk5ytzS7Yyp1Ci0ULcAfpr5F2ezL+CI8j3tfm8zDMz61O5KyZG//jk++e571R9eyKaKiur1hWGU8IxInMuWC27V9STV5WihaAIfTycPTP6X8lQl8Fr6Ltq9dxx+mvWV3rFbJ6/GweOVcvtwwh+zKbdV9KXVzGia4ezC27xQuHT1Nu8tQzYoWihbC4XTy2LQvuOPVcbwfuYG2c27lt9e/YHesVqGs/BgfL/kH3+6aR7bs40CYAxyQioMppHHxkOmcl/Hj0y9IqSZKC0ULEh4ewZPTFnDbaxfwL7OMdh/8hp9f+YjdsVqkvIO5fLT0GVYXLGVD2BGOOh1EuAwDKttyScRwLj93Jqm9MuyOqVSD0ELRwrSNbMdT//U5t86ZyEtHPqPNJ1HMmHS/3bFahFUbv+KzVa+SXbaWzdb1DR3DvGS4YxkefyFTzr9dO2tULZKeHttC7S/cw8y5PyEnzMv5lfHcc9mL9OiaYnesZqWs/BiffP0iK3Z9xiazlz3hvi4zelQYBjp6cl6fK7S9QbUYdZ0eq4WiBTtQtJeH5s7gS9de2nsNV7Ybz6yr/66nz9bBd++GF1lfksWm8KMcczgIM4a0E+EMaDuQiRk/Y1T6xXbHVKrBaaFo5RZ89y7Pr/0zWyK8pJ1w8vNh/8u4zCl2x2oSKipO8MW3b7Bsx0dsqtxZfZZSnNvLQG88QztfwOSxM0mISbQ5qVKhpYVC4XZX8uT7v+KDY19x1CFMcHfn3qteaXVfgF6Ph+82zOer9e+w+Wg2W8OO+e74ZgwpFU76h6dyfpreu0G1PlooVLVd+Vt5+N+3sCS8mAS3l5/GTeaWSX9q0V+KW3et4YvvX2dD8fdsdRZT4HIA0MntpZ8njoHx53DZqFu0y3bVqmmhUKf46KuXeGnr38kNh8Hl4cw672FGDJxgd6wGsXXXOpasfY8NB79lm9nHLquTvSiPlzR3e/q1T2d8+nVkDhjfogukUmdCC4UKqPxEGX9793Y+rszCjTCksh2DO57D5DF3NJs+h8rKj7Fk1Yes3jmfnGNbyXUcYb91Q59wryG1MpzUiBTOTbmcCzOv1jOUlKqFFgpVp005K3lu4W9Y5zhAkeuH4/UDIvoxYdANjM24rEn88vZ6PGzPW8+StXPZUvA9uzz55IS5OeHwFYY4t5dkd3uS2vRmSI9xjB9xDdHtY21OrVTzoIVCBcXtrmTBinf4asu7bHTvIMc6ZNOl0jDAdOGcnj9i8tiZtGsbFdIc5SfKWLd1Gdm5y9hzaCP7y/dSwBH2udwccfraF1zGkFzhoKejM6mxQzlv0BUM6nNOkyhoSjVHWihUvWzYkcW/V7zA+iOr2BRezgmH0NbrpU9FBFHSlvaO9nQI70h0ZAKx7bvSOSaJbgl96NW1H1HtOp60rIqKExws2UfhoTyKD+/n0NGDHCkr4ujxYo5VHObQiQIOVh7koKOMfWGGSpHq18a6vXTxhJMgHekS2YNB3ccwLvMq7XVVqQak96NQ9TKwTyYD+/g+N4cOF/Dh18+Tlb+QfCnmoOMQJc4STpi9cBzfowDY6ntte4+XKK9QKYYyB5Q5HHWuy+kwdHVCJ08bBngS6Nq+N327DmVY2gS6d0oK5dtUSp2G7lGoevN6PBQdOcCu/K3sLdxOweE9HDq2n8PlhZS6D3PMewyXuGjjaEOksy1tXe1pF96BdhEdiWoTR8d2CcR06ERcdCLdO/fRW7kqZSPdo1Ah4XA6SYhJJCEmkUzG2R1HKRUidR8PUEop1eppoVBKKVUnLRRKKaXqpIVCKaVUnbRQKKWUqpMWCqWUUnXSQqGUUqpOWiiUUkrVqcVdmS0iBcCus1hEPFDYQHFCrTllheaVtzllheaVtzllheaV92yy9jLGJASa0OIKxdkSkazaLmNvappTVmheeZtTVmheeZtTVmheeUOVVQ89KaWUqpMWCqWUUnXSQnGqF+wOcAaaU1ZoXnmbU1ZoXnmbU1ZoXnlDklXbKJRSStVJ9yiUUkrVSQuFUkqpOrWaQiEil4jIFhHZLiL3Bpj+axHZKCLrRGShiPTym+YRkTXW4+Mmkne6iBT45brZb9o0EdlmPaY1gayP++XcKiIlftMadduKyCsiclBEsmuZLiLypPVe1onIML9pjbpdg8x7g5VzvYgsE5EhftNyrfFrRCTkt30MIus4ETns9/e+329anZ8hG7L+xi9ntvU5jbWmNep2tdbZQ0QWWd9RG0TkVwHmCd1n1xjT4h+AE9gB9AbCgbXAgBrzjAfaWsO3A2/7TTvaBPNOB54O8NpYIMf6N8YajrEza435fwm8YuO2PR8YBmTXMv3HwGeAAKOA7+zYrmeQd3RVDuDSqrzW81wgvglt23HAp2f7GWqMrDXmnQR8add2tdbZFRhmDUfhuzt9ze+EkH12W8sexUhguzEmxxhTAbwFTPafwRizyBhTZj39FujeyBn9nTZvHX4EzDfGFBtjDgHzgUtClBPOPOv1wJwQ5qmTMWYJUFzHLJOB143Pt0BHEelK42/XoPIaY5ZZecDmz20Q27Y2Z/N5r5czzGrrZxbAGLPPGLPKGi4FNgHdaswWss9uaykU3YA9fs/zOHUj+/tvfJW5SqSIZInItyJyRQjy1RRs3qusXcz3RKTHGb62oQS9PutwXjLwpd/oxt62p1Pb+2ns7VofNT+3BviPiKwUkVttylTTuSKyVkQ+E5GB1rgmu21FpC2+L9X3/Ubbul1FJAkYCnxXY1LIPruuM07ZwonIVCATuMBvdC9jzF4R6Q18KSLrjTE77ElY7RNgjjHmhIjcBrwGXGhzptO5DnjPGOPxG9cUt22zIyLj8RWK8/xGn2dt207AfBHZbP2StssqfH/voyLyY+BDIMXGPMGYBHxjjPHf+7Btu4pIe3xFa5Yx5khjrBNazx7FXqCH3/Pu1riTiMhE4D7gcmPMiarxxpi91r85wGJ81TyUTpvXGFPkl/ElYHiwr21gZ7K+66ixC2/Dtj2d2t5PY2/XoInIYHyfgcnGmKKq8X7b9iAwF98hHtsYY44YY45aw/OAMBGJpwlvW+r+zDbqdhWRMHxF4g1jzAcBZgndZ7cxG2TseuDbc8rBd9ijqrFsYI15huJrUEupMT4GiLCG44FthL6hLZi8Xf2GpwDfmh8arnZauWOs4Vg7s1rzpeFrBBQ7t621riRqb3D9CSc3CK6wY7ueQd6ewHZgdI3x7YAov+FlwCU2Z+1S9ffH9+W629rOQX2GGjOrNT0aXztGuyawXQV4HXiijnlC9tltFYeejDFuEfkF8AW+MyxeMcZsEJE/AlnGmI+BR4D2wLsiArDbGHM50B/4h4h48e2BPWiM2dgE8t4pIpcDbnwf5unWa4tF5E/A99bi/mhO3m22Iyv4fpm9ZaxPrqXRt62IzMF39k28iOQB/wuEWe/leWAevrNHtgNlwE3WtEbdrmeQ934gDnjW+ty6ja/30M7AXGucC3jTGPO5zVmvBm4XETdwHLjO+jwE/AzZnBV8P8D+Y4w55vfSRt+uljHAz4D1IrLGGvc/+H4ohPyzq114KKWUqlNraaNQSilVT1oolFJK1UkLhVJKqTppoVBKKVUnLRRKKaXqpIVCNXkicjSIeWZZ3S001DqvEJEBDbi8ZWfx2qPWv4ki8l4d83UUkTvqux6laqOFQrUUs4AzKhQi4qxj8hVAgxUKY8zoBlhGvjHm6jpm6QhooVANTguFajas+xkstjpB3Cwib1h98N8JJAKLRGSRNe/FIrJcRFaJyLtWHzlV9xJ4SERWAdeIyC0i8r3VUd37ItJWREYDlwOPWPcc6CMiGVbHhetEZK6IxFjLWyy++21kicgmERkhIh9Y/f7/2S/7Ub/he8R3P4O1IvJggPeZbGVfX2MZSWLdP0FEBorICivfOhFJAR4E+ljjHhGR9uK7t8oqa1mT/ZazSUReFN+9Df4jIm2saX1FZIGVbZWI9LHG/8baTutE5IEG/cOqpi/Ul57rQx9n+8C6ZwW+K2kP4+urxgEsx9dBG/jdIwBfdyBLsLpeAO4B7veb77d+y47zG/4z8EtreDZwtd+0dcAF1vAfsbpSwNc/1UPW8K+AfHz3DojA10tnXI33cCm+bh+q7n1ySlcKwMfAjdbwz/1em4TV5QTwFHCDNRwOtKFGlxT4rhzu4LdNtuPr3iEJ3xX9Gda0d4Cp1vB3wBRrOBLfXtrFwAvWax3Ap8D5dn8u9NF4j1bRhYdqUVYYY/IArK4MkoClNeYZhe+w0TdWVwvh+IpKlbf9hgdZv9o74uvC5YuaKxSRaKCjMeYra9RrwLt+s1R1U7Ie2GCM2We9LgdfZ2xFfvNOBF411r1PTOCuFMYAV1nD/wQeCjDPcuA+EekOfGCM2Wa915OiA38RkfMBL76upTtb03YaY9ZYwyuBJBGJAroZY+Za2cqt93ExvmKx2pq/Pb5eX+3siVY1Ii0Uqrk54TfsIfBnWPDdqOX6Wpbh33fPbOAKY8xaEZmOb6+lvpm8NfJ5a8kXjDr71jHGvCki3+HrCG6e+Lqaz6kx2w1AAjDcGFMpIrn49hL8M4NvO7apY3UC/NUY848zyK9aEG2jUC1FKb5bRILvTm9jRKQvgIi0E5HUWl4XBewTXxfONwRanjHmMHBIRMZa034GfEX9zAduqjpDS6z7MNfwDb5OFKmRqZr47t+RY4x5EvgIGMzJ2wB8vZ8etIrEeKDXqUv6gfHdOS1PrBtIiUiElfMLYIZfO0838d2LQbUSWihUS/EC8LmILDLGFODrTXeOiKzDd5gmrZbX/V98x+W/ATb7jX8L+I2IrLYadKfha9xeB2Tga6c4Y8bX0+jHQJZ16OzuALP9Cvi5iKyn9juRXQtkW8sYhO8WmEX4Drdli8gjwBtAprWcG2u8v9r8DF/PxOvwtaV0Mcb8B3gTWG4t6z1OLkiqhdPeY5VSStVJ9yiUUkrVSQuFUkqpOmmhUEopVSctFEoppeqkhUIppVSdtFAopZSqkxYKpZRSdfr/fx29TeTNSRwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "for value, bootstrap in conditions.items():\n",
    "    plt.plot(results[f'points_{bootstrap}'], results[f'energies_{bootstrap}'], label = f'{bootstrap}')\n",
    "plt.plot(results['points_np'], results['energies_np'], label = 'numpy')\n",
    "plt.legend()\n",
    "plt.title('Dissociation profile')\n",
    "plt.xlabel('Interatomic distance')\n",
    "plt.ylabel('Energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare number of evaluations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no bootstrapping\n",
      "Total evaluations for no bootstrapping:\n",
      "2931\n",
      "bootstrapping\n",
      "Total evaluations for bootstrapping:\n",
      "991\n",
      "np\n",
      "Total evaluations for np:\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "for condition, result_full in results_full.items():\n",
    "    print(condition)\n",
    "    print('Total evaluations for ' + condition + ':')\n",
    "    sum = 0\n",
    "    for key in result_full:\n",
    "        if condition is not 'np':\n",
    "            sum += result_full[key]['raw_result']['cost_function_evals']\n",
    "        else:\n",
    "            sum = 0\n",
    "    print(sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extrapolation \n",
    "\n",
    "Here, an extrapolator added that will try to fit each (param,point) set to some specified function and suggest an initial parameter set for the next point (degree of freedom). \n",
    "\n",
    "- Extrapolator is based on an external extrapolator that sets the 'window', that is, the number of previous points to use for extrapolation, while the internal extrapolator proceeds with the actual extrapolation.\n",
    "- In practice, the user sets the window by specifying an integer value to num_bootstrap - which is also the number of previous points to use for bootstrapping. Additionally, the external extrapolator defines the space within how to extrapolate - here PCA, clustering and the standard window approach. \n",
    "\n",
    "In practice, if no extrapolator is defined and bootstrapping is True, then all previous points will be bootstrapped. If an extrapolator list is defined and no points are specified for bootstrapping, then the extrapolation will be done based on all previous points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Window Extrapolator: An extrapolator which wraps another extrapolator, limiting the internal extrapolator's ground truth parameter set to a fixed window size\n",
    "2. PCA Extrapolator: A wrapper extrapolator which reduces the points' dimensionality with PCA, performs extrapolation in the transformed pca space, and untransforms the results before returning.\n",
    "3. Sieve Extrapolator: A wrapper extrapolator which performs an extrapolation, then clusters the extrapolated parameter values into two large and small clusters, and sets the small clusters' parameters to zero.\n",
    "4. Polynomial Extrapolator: An extrapolator based on fitting each parameter to a polynomial function of a user-specified degree.\n",
    "5. Differential Extrapolator: An extrapolator based on treating each param set as a point in space, and performing regression to predict the param set for the next point. A user-specified degree also adds derivatives to the values in the point vectors which serve as features in the training data for the linear regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Here we test two different extrapolation techniques"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define different extrapolators\n",
    "degree = 1\n",
    "extrap_poly = Extrapolator.factory(\"poly\", degree = degree)\n",
    "extrap_diff = Extrapolator.factory(\"diff_model\", degree = degree)\n",
    "extrapolators = {'poly': extrap_poly, 'diff_model': extrap_diff}\n",
    "\n",
    "for key in extrapolators:\n",
    "    extrap_internal = extrapolators[key]\n",
    "    extrap = Extrapolator.factory(\"window\", extrapolator = extrap_internal)\n",
    "    # define extrapolator\n",
    "    # BOPES sampler\n",
    "    bs_extr = BOPESSampler(\n",
    "        gss=me_gsc\n",
    "        ,bootstrap=True\n",
    "        ,num_bootstrap=None\n",
    "        ,extrapolator=extrap)\n",
    "    # execute\n",
    "    res = bs_extr.sample(driver, points)\n",
    "    \n",
    "    results_full[f'{key}']= res.raw_results\n",
    "    results[f'points_{key}'] = res.points\n",
    "    results[f'energies_{key}'] = res.energies\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "poly\n",
      "diff_model\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9sklEQVR4nO3dd3hUZfr/8fc9JYWWhN4JHYRAgFA0KiCKuipYURdWImthi6vuT1dd9+uqa8Fdv+raK0YsWNf6xQYCAqIsTYh0MFQRCISE1CnP7485iUOYJEPI5GSS+3Vdc3nm1M8cxrlzznPOc8QYg1JKKVUZh90BlFJK1W9aKJRSSlVJC4VSSqkqaaFQSilVJS0USimlqqSFQimlVJW0UKioIiLPisj/2LDdT0Vkag2X7SoiR0TEWdu5apuI3CciB0Rkb8XcIrJARK6xO6Oqe6L3Uaj6QkSygXaAF/AB64BZwPPGGL+N0Y6L9TmuMcbMtTvL8RCRrsBGoJsxZl+I6QuA14wxL9Z1NmUvPaJQ9c0FxpjmQDdgBnAb8JK9kRoGEXFVM0tXICdUkVCNmxYKVS8ZYw4bYz4CLgemishAABHJFJH7rOHWIvKJiOSKyEERWSQiDmvabSKyW0TyRWSjiIyzxseKyGMissd6PSYisWXbFZGJIrJaRPJEZKuInGONLz/tIiI9ReQrEcmxTtO8LiKJ1rRXCfzgfmydtvmLiCSLiCn7oRaRjiLykZV5i4hcG7T9u0XkbRGZZWX/QUTSKttP1nr/JCLbrCz/CtoHGSKyREQeFZEc4G4RSbDWvV9EtovI30TEISJnAl8CHa3cmRVzh9j2NBFZLyKHRORzEelWk39rVf9poVD1mjFmGbALOC3E5P9nTWtD4JTVXwEjIn2BPwLDraOTs4Fsa5k7gVFAKjAYGAH8DUBERhA41XUrkAicHrRcMAEeBDoC/YEuwN1W3t8AOwgcGTUzxvwzxPJvWrk7ApcCD4jIGUHTJ1jzJAIfAU+GWEewi4A0YCgwEZgWNG0ksI3A/rkfeAJIAHoAo4GrgKut02TnAnus3BlVbVBEJhLY3xcT2P+LgNnV5FRRSguFigZ7gJYhxnuADgTOqXuMMYtMoNHNB8QCJ4mI2xiTbYzZai0zGbjXGLPPGLMfuAf4jTXtt8BMY8yXxhi/MWa3MWZDxY0aY7ZY85RY63iEwI9utUSkC5AO3GaMKTbGrAZeJPCDXWaxMWaOMcYHvEqgoFXlIWPMQWPMDuAx4MqgaXuMMU8YY7xAKXAFcIcxJt8Ykw38b9DnPx7TgQeNMeutdT8ApOpRRcOkhUJFg07AwRDj/wVsAb6wTr3cDoEfcuAmAn/l7xORN0Wko7VMR2B70Dq2W+MgcGSwlWqISDtrnbtFJA94DWgd5mfpCBw0xuRXyNAp6P3eoOFCIK6a9oWdFdbVsZJprQE3x37+4G2Hqxvwb+u0Xy6Bfx+p4bpUPaeFQtVrIjKcwI/P4orTrL+K/58xpgeB0zV/LmuLMMa8YYw5lcAPmgEeshbbY40r09UaB4Ef1Z5hxHrAWmeKMaYFMIXAj2R5tCqW3QO0FJHmFTLsDmO7lelSYV17gt4HZzlA4Cis4uevybZ3AtcbYxKDXvHGmG9qsC5Vz2mhUPWSiLQQkfMJnKt/zRizNsQ854tILxER4DCBU05+EekrImdYjdTFQBFQdnntbOBvItJGRFoDdxE4IoDA1VVXi8g4q4G3k4j0CxGvOXAEOCwinQi0aQT7mUAbwDGMMTuBb4AHRSRORAYROOX1Wqj5w3SriCRZp7VuBN6qZNs+4G3gfhFpbp0m+nMNt/0scIeIDACwGskvq1l8Vd9poVD1zccikk/gL9Y7CZz/v7qSeXsDcwn8aC8FnjbGzCfQPjGDwF/Qe4G2wB3WMvcBy4E1wFpgpTWurOH8auBRAoVnIUf/9V3mHgINx4eB/wP+U2H6gwSKUa6I3BJi+SuBZAJ/+b8P/P0E77n4EFgBrLbyVHU58Q1AAYEG7sXAG8DM492gMeZ9Akdpb1qn37IINIarBkhvuFMqiomIAXpb7TJKRYQeUSillKqSFgqllFJV0lNPSimlqqRHFEoppapUXSdhUad169YmOTnZ7hhKKRVVVqxYccAY0ybUtAZXKJKTk1m+fLndMZRSKqqIyPbKpumpJ6WUUlXSQqGUUqpKWiiUUkpVqcG1USilTozH42HXrl0UFxfbHUVFQFxcHJ07d8btdoe9jBYKpdRRdu3aRfPmzUlOTibQ36JqKIwx5OTksGvXLrp37x72cnrqSSl1lOLiYlq1aqVFogESEVq1anXcR4taKJRSx9Ai0XDV5N9WC4Vl2w9rmHnt3/j46SfsjqKUUvWKFgqLw+mgyHkGB9bsrX5mpVS9lZGRwbvvvnvC61m9ejVz5sypcp7s7GzeeOONE95WOJ599llmzZpVJ9uqSAuFJbnfQJyeI5iSeLujKKXqgRMtFF6vt1bzTJ8+nauuuqpW1xkuLRRBXN6DGF+i3TGUatSys7Pp378/1157LQMGDGD8+PEUFRUBgR/vUaNGMWjQIC666CIOHToUch1z584lLS2NPn368MknnwCBRvqrr76alJQUhgwZwvz58ysdX1payl133cVbb71Famoqb731FgsXLiQ1NZXU1FSGDBlCfn4+t99+O4sWLSI1NZVHH32UzMxMJkyYwBlnnMG4ceM4cuQI48aNY+jQoaSkpPDhhx+Wf8Z+/foxefJk+vfvz6WXXkphYSEQ6IboL3/5CykpKYwYMYItWwLPpLr77rt5+OGHARgzZgy33XYbI0aMoE+fPixatAiAwsJCJk2axEknncRFF13EyJEja6VLI708NoiYHPzOjnbHUKreuOfjH1i3J69W13lSxxb8/YIBVc6zefNmZs+ezQsvvMCkSZN47733mDJlCldddRVPPPEEo0eP5q677uKee+7hscceO2b57Oxsli1bxtatWxk7dixbtmzhqaeeQkRYu3YtGzZsYPz48WzatKnS8ffeey/Lly/nySefBOCCCy7gqaeeIj09nSNHjhAXF8eMGTN4+OGHy4tRZmYmK1euZM2aNbRs2RKv18v7779PixYtOHDgAKNGjWLChAkAbNy4kZdeeon09HSmTZvG008/zS23BJ6cm5CQwNq1a5k1axY33XRT+fqDeb1eli1bxpw5c7jnnnuYO3cuTz/9NElJSaxbt46srCxSU1NP4F/qF3pEEUTceXjcLSkqKLA7ilKNWvfu3ct/5IYNG0Z2djaHDx8mNzeX0aNHAzB16lS+/vrrkMtPmjQJh8NB79696dGjBxs2bGDx4sVMmTIFgH79+tGtWzc2bdpU6fiK0tPT+fOf/8zjjz9Obm4uLlfov7PPOussWrZsCQTuW/jrX//KoEGDOPPMM9m9ezc///wzAF26dCE9PR2AKVOmsHjx4vJ1XHnlleX/Xbp0acjtXHzxxUftH4DFixdzxRVXADBw4EAGDRoUctnjpUcUQRxNSzGlbrK+Wcjws35ldxylbFfdX/6REhsbWz7sdDrLTz2Fq+IloLVxue/tt9/Oeeedx5w5c0hPT+fzzz8POV/Tpk3Lh19//XX279/PihUrcLvdJCcnl9/DUFXGyoaDle0jp9NZ6+0hFekRRZAm7QL/wDu+z7I5iVKqooSEBJKSksrPx7/66qvlRxcVvfPOO/j9frZu3cq2bdvo27cvp512Gq+//joAmzZtYseOHVWOb968Ofn5+eXr3Lp1KykpKdx2220MHz6cDRs2HDNPRYcPH6Zt27a43W7mz5/P9u2/9OS9Y8eO8qOFN954g1NPPbV82ltvvVX+35NPPjnsfZSens7bb78NwLp161i7dm3Yy1ZFjyiCdOjfh/074cjuXLujKKVCeOWVV5g+fTqFhYX06NGDl19+OeR8Xbt2ZcSIEeTl5fHss88SFxfH73//e373u9+RkpKCy+UiMzOT2NjYSsePHTuWGTNmkJqayh133MHixYuZP38+DoeDAQMGcO655+JwOHA6nQwePJiMjAySkpKOyjF58mQuuOACUlJSSEtLo1+/fuXT+vbty1NPPcW0adM46aST+N3vflc+7dChQwwaNIjY2Fhmz54d9v75/e9/z9SpUznppJPo168fAwYMICEh4Tj38rEa3DOz09LSTE1b+Q/nHOC1v64m3r+QaS/8o5aTKRUd1q9fT//+/e2O0aBlZ2dz/vnnk5V17NmLsoevtW7d+rjX6/P58Hg8xMXFsXXrVs4880w2btxITEzMUfOF+jcWkRXGmLRQ69UjiiAJrVrjLs0Ff3O7oyil1HErLCxk7NixeDwejDE8/fTTxxSJmtBCUYHTn4MhqfoZlVKqhpKTk0MeTQDlVzDVRPPmzSPyKGhtzD7GQXzOVnaHUEqpekMLRQUSU4AnJpH9e3baHUUppeoFLRQVuBMCjftrF863OYlSStUPWigqaN45cNpp36Zse4MopVQ9oYWigh7DhwBQvF+fF6yUXbKzsxk4cOAJr+exxx4r72yvMpmZmezZs+eEtxWOU045pU62U9u0UFRw0ohTcfhK8RXGVj+zUqpeO9FC4fP5ajXPN998U6vrqytaKCpwud24PTnga2F3FKUaNa/Xe0w33PPmzWPIkCGkpKQwbdo0SkpKAEKOf/zxx9mzZw9jx45l7Nix+Hw+MjIyGDhwICkpKTz66KO8++67LF++nMmTJ5OamkpRURHJycncdtttDB06lHfeeYcXXniB4cOHM3jwYC655JLywpORkcH06dOP6c48MzOTiRMnMmbMGHr37s0999xT/pmaNWsGwIIFCxgzZgyXXnppeXfjZTc/z5kzh379+jFs2DD+9Kc/cf7559flbg/JlvsoRKQl8BaQDGQDk4wxITuWF5EWwDrgA2PMH+sknz8HI3qJrFJ8ejvsrZ3+gsq1T4FzZ1Q7W8VuuB955BGee+455s2bR58+fbjqqqt45plnmD59OhkZGceMv+mmm3jkkUeYP38+rVu3ZsWKFezevbv8/oXc3FwSExN58sknefjhh0lL++Wm5FatWrFy5UoAcnJyuPbaawH429/+xksvvcQNN9wAhO7OHGDZsmVkZWXRpEkThg8fznnnnXfU+gFWrVrFDz/8QMeOHUlPT2fJkiWkpaVx/fXX8/XXX9O9e/fyXmTtZtcRxe3APGNMb2Ce9b4y/wBC9yUcKc7DeN2t8Ho8dbpZpdQvKnbDPW/ePLp3706fPn2AX7oZ37hxY8jxFfXo0YNt27Zxww038Nlnn9GiReVnDS6//PLy4aysLE477TRSUlJ4/fXX+eGHH8qnherOHAJdjbdq1Yr4+Hguvvjio7oQLzNixAg6d+6Mw+EgNTWV7OxsNmzYQI8ePejevTtAvSkUdt2ZPREYYw2/AiwAbqs4k4gMA9oBnwEh+yCJBEd8ET5/PNvWrqLP0BF1tVml6p8w/vKPlIrdaycmJpKTk1Pj9SUlJfH999/z+eef8+yzz/L2228zc+bMkPMGdxWekZHBBx98wODBg8nMzGTBggWVZix7H0435xW7Uo90V+Enwq4jinbGmJ+s4b0EisFRRMQB/C9wS3UrE5HrRGS5iCzfv3//CYeLbekGYNOyZSe8LqVUzVTshjstLY3s7Ozy0ztl3Yz37ds35HjgqG7ADxw4gN/v55JLLuG+++4rP7VUXVfh+fn5dOjQAY/HU94deZlQ3ZkDfPnllxw8eJCioiI++OCD8iOj6vTt25dt27aVd+NR1t243SJ2RCEic4H2ISbdGfzGGGNEJFQXtr8H5hhjdlX30BFjzPPA8xDoPbZmiX/RqmdnDh2A3B9/PtFVKaVqqGI33I8//jijRo3isssuw+v1Mnz4cKZPn05sbCwvv/zyMeMBrrvuOs455xw6duzIY489xtVXX43f7wfgwQcfBH5plI6Pjw/5NLl//OMfjBw5kjZt2jBy5Mijikqo7swhcFrpkksuYdeuXUyZMuWY9onKxMfH8/TTT3POOefQtGlThg8ffkL7sLbY0s24iGwExhhjfhKRDsACY0zfCvO8DpwG+IFmQAzwtDGmqvaME+pmvMzOzRv56H93E888pj17/wmtS6loo92MhycjI4Pzzz+fSy+99KjxmZmZRz1r+3gdOXKEZs2aYYzhD3/4A7179+bmm2+ujcjljrebcbtOPX0ETLWGpwIfVpzBGDPZGNPVGJNM4PTTrOqKRG3p0rsvLk8+prhJXWxOKaXKvfDCC6SmpjJgwAAOHz7M9ddfb3ck2xqzZwBvi8hvge3AJAARSQOmG2OusSlXOac3B0yi3TGUUvVUZmZmyPEZGRlkZGTUeL0333xzrR9BnChbCoUxJgcYF2L8cuCYImGMyQQyIx4siHAQv7NLXW5SKaXqJb0zuzKufDzulhQVFNidRCmlbKWFohKuZh6Mw8n3X8+zO4pSStlKC0UlmrQL9Mmya+06m5MopZS9tFBUouOAwNW6BT/l2ZxEKaXspYWiEoNGjwPjx5dv14VhSilVP2ihqETzxERiSg9iPM3sjqJUo5OdnU3//v259tprGTBgAOPHj6eoqIgxY8ZQdkPtgQMHSE5OBgKXql544YWcddZZJCcn8+STT/LII48wZMgQRo0axcGDBwEYM2YMN954I6mpqQwcOJBly5bh9/vp3bs3Zd3/+P1+evXqRW10B9RQ6J/LVXD4DwIt7Y6hlG0eWvYQGw5uqNV19mvZj9tGHNMH6DE2b97M7NmzeeGFF5g0aRLvvfdelfNnZWWxatUqiouL6dWrFw899BCrVq3i5ptvZtasWdx0000AFBYWsnr1ar7++mumTZtGVlYWU6ZM4fXXX+emm25i7ty5DB48mDZt2tTGx20Q9IiiKnIIr1OfS6GUHbp3705qaioAw4YNK+8orzJjx46lefPmtGnThoSEBC644AIAUlJSjlq2rOvu008/nby8PHJzc5k2bRqzZs0CYObMmVx99dW1/nmimR5RVEFijuB1JLBv53badulmdxyl6lw4f/lHSsVuuIuKinC5XOWd+hUXF1c6v8PhKH/vcDiO6sI7VBfgXbp0oV27dnz11VcsW7bsmF5iGzs9oqiCOzHwhVr79XybkyilAJKTk1mxYgUA7777bo3WUdZ19+LFi0lISCAhIQGAa665hilTpnDZZZfhdDprJ3ADoYWiCi26BE477du03eYkSimAW265hWeeeYYhQ4Zw4MCBGq0jLi6OIUOGMH36dF566aXy8RMmTODIkSN62ikEW7oZj6Ta6Ga8TNbSxSx8pZQmrrlc/eQDtbJOpeq7htzN+JgxY455PnaZ5cuXc/PNN7No0SIbktWtaOlmPCr0SxuJw1eCKYytfmalVNSaMWMGl1xySfnDjNTRtDG7Ci63G7cnB+NPsDuKUqoWBD/vOtjtt9/O7bfXyeNuopIeUVRD/DkYh14iq5RqvLRQVEOch/G4W+H1eOyOopRSttBCUQ2JL8bvjGPLmpV2R1FKKVtooahGXOsYADZ/+63NSZRSyh5aKKrRulfgcaiHd2gHYUrVV8GdBarap4WiGgNOOx0AT27Dut9EKaXCpYWiGh2798blycOUNLE7ilKNRnZ2Nv369WPy5Mn079+fSy+9lMLCQubNm8eQIUNISUlh2rRplJSUHLXczJkzy3uJBXjhhRe4+eab6zh9w6P3UYTB6c3BGO1uXDU+ex94gJL1tdvNeGz/frT/61+rnW/jxo289NJLpKenM23aNB555BGee+455s2bR58+fbjqqqt45plnjioMkyZN4v777+df//oXbrebl19+meeee65W8zdGekQRBuEgfqcWCqXqUpcuXUhPTwdgypQpzJs3j+7du9OnTx8Apk6dytdff33UMs2aNeOMM87gk08+YcOGDXg8HlJSUuo8e0OjRxRhEHc+pc6WFOQdpmkLvUtbNR7h/OUfKRW7A09MTCQnJ6fa5a655hoeeOAB+vXrpx381RI9ogiDs5kXxMmaRdrduFJ1ZceOHSxduhSAN954g7S0NLKzs9myZQsAr776KqNHjz5muZEjR7Jz507eeOON8ocUqROjhSIMTdo3B2DX2vU2J1Gq8ejbty9PPfUU/fv359ChQ9x88828/PLLXHbZZaSkpOBwOJg+fXrIZSdNmkR6ejpJSUl1nLph0lNPYeic0p992VC4N9/uKEo1Gi6Xi9dee+2ocePGjWPVqlXHzFuxs7/Fixfr1U61SI8owjDotLFgfPiOaF1Vqj7Lzc2lT58+xMfHM27cOLvjNBj6yxeGpi0SiCk9iPE3tzuKUo1CcnIyWVlZx71cYmIimzZtikCixk2PKMLk8B3EoJfIKqUaHy0UYRLHQXwufS6FUqrx0UIRJoktxOtuwd7tW+2OopRSdUoLRZjciYGbf9YuXGBvEKWUqmNaKMKU0LUNAAe27rI5iVJK1S0tFGHqNTwNgOL9JdXMqZSqbXfffTcPP/wwd911F3PnzgVg0aJFDBgwgNTUVIqKirj11lsZMGAAt956a8TzZGZm8sc//vGE54kWenlsmHoPGcEC35eYkji7oyjVaN17773lw6+//jp33HEHU6ZMAeD555/n4MGDOJ1Ou+I1WFoowuRyu3F7cjB+7RRQNR6L3t7EgZ1HanWdrbs047RJfaqd7/777+eVV16hbdu2dOnShWHDhpGRkcH5559Pbm4ub7/9Np9//jmffvop+fn5HDlyhGHDhnHHHXdw+eWXH7O+jIwM4uPjWbVqFfv27WPmzJnMmjWLpUuXMnLkSDIzMwGYPXs2DzzwAMYYzjvvPB566CEAXn75ZR588EESExMZPHgwsbGxAOzfv5/p06ezY8cOAB577LHyXm8bCi0Ux0H8ORiHXiKrVKStWLGCN998k9WrV+P1ehk6dCjDhg0rn37NNdewePFizj//fC699FIg0MX46tWrq1zvoUOHWLp0KR999BETJkxgyZIlvPjiiwwfPpzVq1fTtm1bbrvtNlasWEFSUhLjx4/ngw8+YOTIkfz9739nxYoVJCQkMHbsWIYMGQLAjTfeyM0338ypp57Kjh07OPvss1m/vmH1C6eF4jiI6zClzr54PR5cbrfdcZSKuHD+8o+ERYsWcdFFF9GkSeDJkhMmTKiV9V5wwQWICCkpKbRr1678WRUDBgwgOzub7du3M2bMGNq0CVy8Mnny5PJnXgSPv/zyy8vvAJ87dy7r1q0r30ZeXh5HjtTuUZjdbCkUItISeAtIBrKBScaYQyHm6wq8CHQBDPArY0x2nQWtmKdJCX5vLBuWf8fAk0+1K4ZSqobKThc5HI7y4bL3Xq8Xdw3+APT7/Xz77bfExTXc9ku7rnq6HZhnjOkNzLPehzIL+Jcxpj8wAthXR/lCimsV+GJtXf5fO2Mo1eCdfvrpfPDBBxQVFZGfn8/HH39cJ9sdMWIECxcu5MCBA/h8PmbPns3o0aMZOXIkCxcuJCcnB4/HwzvvvFO+zPjx43niiSfK31d3+isa2XXqaSIwxhp+BVgA3BY8g4icBLiMMV8CGGNsP5Zr26cbB3+GvJ3VP2VLKVVzQ4cO5fLLL2fw4MG0bduW4cOH18l2O3TowIwZMxg7dmx5Y/bEiROBwCW6J598MomJiaSmppYv8/jjj/OHP/yBQYMG4fV6Of3003n22WfrJG9dEWNM3W9UJNcYk2gNC3Co7H3QPBcC1wClQHdgLnC7McYXYn3XAdcBdO3addj27dsjknvfzu28c/9W4v3zmPb8/RHZhlJ2W79+Pf3797c7hoqgUP/GIrLCGJMWav6InXoSkbkikhXiNTF4PhOoVKGqlQs4DbgFGA70ADJCbcsY87wxJs0Yk1bW2BQJbbt0w1V6GFPaLGLbUEqp+iZip56MMWdWNk1EfhaRDsaYn0SkA6HbHnYBq40x26xlPgBGAS9FIm+4XL4cQB+vqFR9df/99x/VhgBw2WWXceedd9qUKPrZ1UbxETAVmGH998MQ8/wXSBSRNsaY/cAZwPK6i1iZg/gdyXaHUCqijDEEzgpHnzvvvFOLQhVq0txg11VPM4CzRGQzcKb1HhFJE5EXAay2iFuAeSKyFhDgBZvylhP3EUpjksjPzbU7ilIRERcXR05OTo1+UFT9ZowhJyfnuC/lteWIwhiTAxzzQFtjzHICDdhl778EBtVhtGo5m3uhyMmahfNIn3iJ3XGUqnWdO3dm165d7N+/3+4oKgLi4uLo3LnzcS2jd2Yfp6YdWnBkG+xZtzFwka9SDYzb7aZ79+52x1D1iHYzfpy6pg4EoHCv7bd1KKVUndBCcZxS0scifi++I3owppRqHPTX7jjFN22K23MQ429hdxSllKoTekRRAw5fDoaWdsdQSqk6oYWiBsR5CK+7LV6Px+4oSikVcVooasCZUIjP1YRvPn7X7ihKKRVxWihqoFNa4GEu2d9k2ZxEKaUiTwtFDZx+6RU4vUV4D8TbHUUppSIurEIhIv8rIgMiHSZaxMTF4S7djt8c392NSikVjcI9olgPPC8i34nIdBFJiGSoaCAxeymJ68Te7VvtjqKUUhEVVqEwxrxojEkHriLwnOs1IvKGiIyNZLj6rElnJ4iTJe+8Z3cUpZSKqLDbKETECfSzXgeA74E/i8ibEcpWr6VecA4AeVsO25xEKaUiK9w2ikeBjcCvgAeMMcOMMQ8ZYy4AhkQyYH3Vb9hIYkp+xl/U2u4oSikVUeF24bEG+JsxpiDEtBG1mCeqOPw78Dr74PV4cLnddsdRSqmICPfU0/dAXxEZGvTqKSIuY0yjPffibHEYb0wCy+d+ancUpZSKmHALxdPAt8DzBJ4ytxR4B9goIuMjlK3ea5MSuDx2y4JlNidRSqnICbdQ7AGGGGPSjDHDCLRLbAPOAv4ZqXD13emTfo3DV0rpfj3tpJRquMItFH2MMT+UvTHGrAP6GWO2RSZWdGiemEhM6Q78Pr3xTinVcIVbKNaJyDMiMtp6PW2NiwUadReq4txNaWxnDu3fa3cUpZSKiHALxVRgC3CT9doGZBAoEo32pjuA2A5gHG6+fnO23VGUUioiqr081rrRbo4xZizwvyFmadQPjx5w1mkseRMOrt9vdxSllIqIao8ojDE+wK/9O4WWOuZM3CUH8RfoE++UUg1TuDfcHQHWisiXQPlNd8aYP0UkVZRx+rfjc3SzO4ZSSkVEuIXiP9ZLheBomkOxfwjfL/qKwaedYXccpZSqVWEVCmPMKyISD3Q1xmyMcKaok9S3DYXr4YfPF2qhUEo1OOF2CngBsBr4zHqfKiIfRTBXVDnt8kmI30vxT3YnUUqp2hfu5bF3E+j8LxfAGLMa6BGRRFGoVftOxJTsxHg62h1FKaVqXbiFwhOi8z9/bYeJZg7Hbkpju1GQ12j7SFRKNVDhFoofROTXgFNEeovIE8A3EcwVddxtS/E7Y1j41ht2R1FKqVoVbqG4ARgAlACzgTwCd2grS6/T0wDYt2anzUmUUqp2hfvM7EJjzJ3GmOFWD7J3GmOKIx0umgwffx6u0jx8eXpfolKqYQnr8lgR6QPcAiQHL2OM0WtBLS63G5c3G7+zq91RlFKqVoV7w907wLPAi4AvcnGimzTZTymD2LRyGX2GNtonxCqlGphwC4XXGPNMRJM0AAk9EyjaCis//kwLhVKqwQi3MftjEfm9iHQQkZZlr4gmi0Lpky4B46dwZ6N+RIdSqoEJ94hiqvXfW4PGGfSmu6O079aT2OLFGNPe7ihKKVVrwu3rqXukgzQUDsdOSp1DKC0uJiYuzu44Sil1wqo89SQifwkavqzCtAciFSqauZIK8bniWfTem3ZHUUqpWlFdG8UVQcN3VJh2Tk03arVxfCkim63/JlUy3z9F5AcRWS8ij4uI1HSbdaXrqAEA7PrvJpuTKKVU7aiuUEglw6HeH4/bgXnGmN7APOv90SsXOQVIBwYBA4HhwOgT2GadOGXCJTi9Bfhym9kdRSmlakV1hcJUMhzq/fGYCLxiDb8CXFjJtuOAGCAWcAM/n8A260RMXBzu0u346WJ3FKWUqhXVFYrBIpInIvnAIGu47H3KCWy3nTGm7OkNe4F2FWcwxiwF5gM/Wa/PjTHrQ61MRK4TkeUisnz//v0nEKt2SOxeSmI7sGPTOrujKKXUCauyUBhjnMaYFsaY5sYYlzVc9t5d1bIiMldEskK8JlbYhiHE0YmI9AL6A52BTsAZInJaJTmft/qgSmvTpk01HznymiXHgTj47r0P7I6ilFInLNz7KI6bMebMyqaJyM8i0sEY85OIdAD2hZjtIuBbY8wRa5lPgZOBRREJXIvSLpzAp08cIP/HArujKKXUCQv3zuza9hG/3MQ3FfgwxDw7gNEi4hIRN4GG7JCnnuqbHgMGEVP8E6b4mDNqSikVdewqFDOAs0RkM3Cm9R4RSRORF6153gW2AmuB74HvjTEf2xG2JhxmB153Ml6PduehlIpuETv1VBVjTA4wLsT45cA11rAPuL6Oo9UaZ2I+xSXNWPrx+5x28SS74yilVI3ZdUTR4HVIDfR68uM339ucRCmlTowWigg5fdKvcfiK8eRof09Kqehmy6mnxiC+aVNiSrZjRG+8U0pFNz2iiCCJ+YnS2E7s27nd7ihKKVVjWigiKL6TE+Nwsvidt+2OopRSNaaFIoKGTDgHjJ/cdXl2R1FKqRrTQhFB/YaNJK5oMz7fQL2fQikVtbRQRJir1U5KY9sw59mn7I6ilFI1ooUiwkZf/2vE7+HASj39pJSKTlooIiy530BiS7LwOAZRkHfY7jhKKXXctFDUgfgueXjdLfjkscftjqKUUsdNC0UdOPfG3+H0FlKwNdbuKEopddy0UNSBpDbtifF+T0lMCnu3b7U7jlJKHRctFHUkYYAbvzOWL5962e4oSil1XLRQ1JELbvgT7tKDlO5ta3cUpZQ6Lloo6khMXBwuWUtJXH82rPjO7jhKKRU2LRR1qNOpnTEOJ99mhnryq1JK1U9aKOrQ2VdfS2zxbnx53e2OopRSYdNCUccccRsoju/J0o/ftzuKUkqFRQtFHTtp4kgA1n+83OYkSikVHi0UdWzUuROIK9qCr/Qk7VFWKRUVtFDYwJnwI6VxHfgi80W7oyilVLW0UNjglIyLEL+XvUt/tjuKUkpVSwuFDfoMHUFs8Xq8DKKooMDuOEopVSUtFDaJ7bAfT0winzyhPcoqpeo3LRQ2Gf/Ha3H4islbb3cSpZSqmhYKm7Tt0o2Y0rV43IPI2bvb7jhKKVUpLRQ2at7Hg88Vz2ePP2t3FKWUqpQWChudd8MNuErzKN7d0u4oSqkoV1RQwJ4fN0dk3a6IrFWFpWmLBNxmDcWxI9mydhW9UobYHUkpVQ8UFRSwdc0qftq8kcO7f6b4cBHeAh/+Egd43RhfHJgmGGmCcTTF52yKzxlPXPGP/PaV3rWeRwuFzdoMT2LHGjdLXnyHXv/WQqFUQ1VaXMzWNcvJXpPF4d37KDlUjL9IMKWxGH88mGYYRzN8zuZ4Xc1AHEB362VxgMNZjNMU4PAXIKYAh/8QDoqIkVJcbfwRya6FwmbnXjudl697F29hV7ujKKVq6ND+vfywZBF7N26hcP8RvEeAkliMrxlIC3yORLzuFhiHC+hjvQABh6sYlzcfhz8f8efgYgduKUbivLibu4hLakpC+za07d6drv1OIqlN+zr/fFoobOZyu3G61lEUM5oV8z5l2Lhz7Y6klKoge0MWG5cu5WD2bkpySvEXujG+5hiS8DuT8LhbgLQCWpUv43QW4TSHcfgP4/RvxeUrRGI8xLRw0bRdEu16dqfn4KG07VL//0jUQlEP9D43hTVfOFjz7iItFErZwOvxsO6/S9j23QrydubgzXNhPM0xtMLrboPP1RToab1AXB7c5hAO3yGcvo24nIU4m/lo0q45bXt1o/ew4XTsXvttBXbRQlEPnHbxJDZ+OBOfnERpcTExcXF2R1KqQcrekMX3X37J4W0HrGKQiHG0wuNug98ZC1jthE4/Mb6DOHz7ifF9j7gLcbdw0qJzK7qlDqL/8NGN6v9TLRT1REybreQfGcebf/sHVz18v91xlIpaXo+HrCUL2bx0GUd25+MvaAr+1nid7fDGJACDARCnB7cvB4f/ALHebCS2iPi28XQc2JvBo88ioVVrez9IPaKFop644r7/Ydb1b1BcOoy927fSvltPuyMpVe/t2LSOlf/3KbnbcvAfaYLxt8fj7oDPFQ+MAgJtBW7/Xty+jbjJI65tDJ0G9SFt/Hk0bZFg7weIEloo6omYuDhaD8ll94Zk5tyXybQX/mF3JKXqjdLiYr6d8z47l6+jZD+Yklb4nR0pjW1N2ekip7MQt283bu9K4uKKaNqpOT1HpTHo1PG43G57P0CUE2OM3RlqVVpamlm+PHofMzpz6sOUxAxk2GVeRow/3+44StW50uJivvnoPXYuW4fnYAzG14HSmC5WGwJg/MSU7MNh9iCxh4jvEEOv9DSGjNWCcCJEZIUxJi3UNFuOKETkMuBuoD8wwhgT8pddRM4B/g04gReNMTPqLKRNBk7ux4p3YO3r67VQqAbP6/GwbM5HbPtmNaX7HRhfBzzuLvhcHYAOOJylxPh2Euv9L85mxbTq25GREybSpuOZdkdvVGw5ohCR/oAfeA64JVShEBEnsAk4C9gF/Be40hizrqp1R/sRBcDMa/6HItdYOvVbzYU3/dnuOErVmh2b1vHtO+9zJLsIU9IeT0w369JTEL+HmJLdOGQPzqRCOqT2IP2iy7QdoY7UuyMKY8x6ABGparYRwBZjzDZr3jeBiUCVhaIhOOfOq/jkH2s5sCpRL5dVUcvr8bD04//w45I1eHKaYOhCSWxHkJPB4SfW8RMx3jU4mxbQZmAX0i++jIRWZ9sdW4VQnxuzOwE7g97vAkaGmlFErgOuA+jatf7f5Vidjt17E5eUSX6BXi6rokd+bi7zXnmZg+sPYArb4nUn43W3AcbhdBXhLv2ReP9mmiXHMeqyi+jaR08fRYuIFQoRmQuE6pTkTmPMh7W5LWPM88DzEDj1VJvrtssV9+vlsqp+O7R/Lwtee41D63MxxR0pje2O3xm4RyHG8TNu33pim+XRaVhPRk+arEfGUSxihcIYc6J/LuwGugS972yNaxT0cllV3+Ts3c38Wa+St7kAU9qJ0phk/M6h4IBYx25ivMuJa13K4AlnM/DkK+2Oq2pRfT719F+gt4h0J1AgrgB+bW+kunXhTX+2LpdNZ9kXn+hVUKpOFRUUMDfzRfZ/vw9T3IWSuO4Yx4hA+wK7ifV9R1w7P2kXn0efob+xO66KILsuj70IeAJoA/yfiKw2xpwtIh0JXAb7K2OMV0T+CHxO4PLYmcaYH+zIaye9XFbVFa/Hwzcfv8u2+Wvx5bXB4+6Nz5USOGKQncR6v6VJFxh+6QR6pWj7QmOiN9xFAb1cVkXK+mVLWP7ep5TubYrP0RtPbOCxvO6SgzjNZmLbFTDsknPpPyLd5qQq0urd5bHq+Ojlsqq2FOQd5suZL3EwKxe/pwcl8V2BM3C6C3GXbibGvZKeY1M5+YKL9C5nVU4LRRTQy2XViVizZD6r3v0S78FWeNx98LlSweknzvMj8f6vaDO4DWdmXEN8Uz21qULTQhEl9HJZFa6CvMN8/uLz5P6Qj9/bk5L4LsCZuJ25uD1ZNGtVyKgpF2s7gwqbFooooZfLqqps+X4537z2Pt59SZS6++FzDQOXjzjPNpqYebQf3oWzpv4Wl/tiu6OqKKSN2VGmrHfZIReXMurcCXbHUTbxejx8/e5stn+9AX9BN4rje4A4cXnycHnXE9epkPSrLie530C7o6oooY3ZDUjgclk/WW/l0bLdMvoMHWF3JFVHDucc4Ivnnid/sxcffSmN7Qx0JlZ2Eu/7mtaDW3HOtdOJibvQ7qiqgdFCEWVGjD+fHaseZN/WoXz97yzi/5ZAl9597Y6lImTL2lV8M+tdPPta4nH3w+cahcNVSkzJJpq6VzNgwmkMP2uq3TFVA6ennqLU7Lvu5uDP6cQVbWXCjPNo07FL9QupqLD0/z5gw/8tw3+kCyVxvTAOJ67SPNz+dcR1LeWMa66ifdcedsdUDUxVp560UESxV2+5k7z8scQV/cCkx6fSPDHR7kiqBrweD5+99Bz7lu/D7+1FSVxnAGKLf8Lh3Eir1CTOvfZ3ev+MiigtFA1Y5g1/pcBzJnFFK/jNczfoj0mUOJxzgM+eeZaCLeB19McTkwTGT1zRVhxNt9PzrEGcfvEVdsdUjYgWigbu5evvpFDGEV/8DVe9cJveUVtPZW/IYtHMNyn9ORGPuz8+VzwOXwkxpRtwt9rPiMkT6Tcs5CNXlIo4LRQNnNfjYdb0+yhyjybeu4BpL95rdyRlWTHvU9a+vwhfXidK4npjHC5cnjzcvnXEd/Nw5vXTtH1J1QtaKBoBr8fDK9f+k+K4k2nCPK5+Vrv5sIPX42Hea5nsXrodf0kvqy8liCn+GadjAy0HJXD2tdcT37SpzUmVOpoWikaiqKCAN37/NMXxw2gaO5eMfz9gd6RGIT83l8+eeZr8TV589KM0tjUAsUU/4ojfRvKYfpxxhT6vQdVvWigakfzcXN7+0yyK40+iReICfvPQfXZHapC2/bCGJbPepnRvgtXe0ASHr5SYko24kn5myKSzGJQ+1u6YSoVNC0Ujs3/PTj66Yw4lsT1I6vANV97zd7sjNQhLP36fDZ//F39e0P0Nnnxc3nXEdy0O3N+gnTWqKKWFohHasWkdX9z/HSVxXYgr+Y7RN59Nr8EhvwOqEvm5uXz+3DPkbSzB5+9NaVwHAGKKf8Lp3ETLFG1vUA2HFopGaufmjXzx4JsUu0/B5S0krsV3TP7nPXr5bBXKn92Q0xpPTF98rnjE7yG2eAuOZrvpeWYKp198pd0xlap1Wigauc9mPs/uhUJxfE/iirbQdZyTs676rd2x6oWiggK+nPkCB9bk4C/tWX6Vkrs0F5d/PXGdSxnz2yl07N7b5qRKRZYWCkVpcTGz77iX4rxReF3xxHuWcM7/TG50P4Bej4clH77Djwuz8OW1wePujc/VJHBXdPGPSGw2HUaUPbtBj7xU46GFQpXbtHIZix//kqK4k3GXHqJp+9Vcfs9dDfpHcc2S+az+YC6efc3wOXvjiWkJgLvkIE6zmZh2BYy8fIJ22a4aNS0U6hgf/PtRDqxsRUl8Z+IK19H3knacOvEyu2PVijXfLCRrzlcU73Hg9/egJK4TAE5vAW7PJpwtDtBz7BBOvuCiBl0glToeWihUSAV5h3n7jocoKUnH73ARW5KFu/VBRvx6QtT0OZSfm8vXb73OvrW78Oe3xO/oVn7DW6AReivSZDfthnTkjCkZeoWSUpXQQqGqtHrBXFZmfovXcRKemMRfztfHZ9P11N6MuXxyvfjL2+vxsG7ZYrI+W0DxHjCejpTGdsPvjAECDdBO73YcTfeT1KcVp115Ba3ad7I5tVLRQQuFCktpcTFfvvISPy/fe9QVQDElB3CygWa9nJw9/XoSWrWOaI6CvMOsnPs5u75fT/G+EvyFzYHWeFzt8LmbASB+LzElu3A4dhHTtpQ+40YxdOzZ9aKgKRWNtFCoGln51Res+XABvkPtKY3tg98Zg9NbhLs0G+QI4ixCYktxxkNMYhzN27emfc/u9Bg0hKQ27Y9aV1FBAT9lb2XfjmwO7/2ZgkO5lBwuxFtQiq/E4C9yYUoTMY62lMa0wTh++cF3lR7G5f0ZcR5A4gtJ6t2K9Msnaa+rStUiLRTqhO3buZ15L2ZSuN2FMR3wO5rjczUvP+1TkdNbiNNXgBEXPmc8fmc1D1QyPmJK9uPw70Pch3C28JDYow0pZ4yjx4BBEfhESqlgVRUKV12HUdGpbZdux/QZ5fV42Lcrmx/XrGF/9g4KD+TiyfPgK3KALwbjjwe8uCgFvIjbhyMWXPFOYprFE5fYnBZtW5HUsRPdTxqqj3JVqp7SQqFqzOV207F770Z3055SjY3D7gBKKaXqNy0USimlqqSFQimlVJW0UCillKqSFgqllFJV0kKhlFKqSloolFJKVUkLhVJKqSo1uC48RGQ/sP0EVtEaOFBLcSItmrJCdOWNpqwQXXmjKStEV94TydrNGNMm1IQGVyhOlIgsr6y/k/ommrJCdOWNpqwQXXmjKStEV95IZdVTT0oppaqkhUIppVSVtFAc63m7AxyHaMoK0ZU3mrJCdOWNpqwQXXkjklXbKJRSSlVJjyiUUkpVSQuFUkqpKjWaQiEi54jIRhHZIiK3h5j+ZxFZJyJrRGSeiHQLmuYTkdXW66N6kjdDRPYH5bomaNpUEdlsvabWg6yPBuXcJCK5QdPqdN+KyEwR2SciWZVMFxF53Posa0RkaNC0Ot2vYeadbOVcKyLfiMjgoGnZ1vjVIhLx5wOHkXWMiBwO+ve+K2hald8hG7LeGpQzy/qetrSm1el+tbbZRUTmW79RP4jIjSHmidx31xjT4F+AE9gK9ABigO+BkyrMMxZoYg3/DngraNqRepg3A3gyxLItgW3Wf5Os4SQ7s1aY/wZgpo379nRgKJBVyfRfAZ8CAowCvrNjvx5H3lPKcgDnluW13mcDrevRvh0DfHKi36G6yFph3guAr+zar9Y2OwBDreHmwKYQvwkR++42liOKEcAWY8w2Y0wp8CYwMXgGY8x8Y0yh9fZboHMdZwxWbd4qnA18aYw5aIw5BHwJnBOhnHD8Wa8EZkcwT5WMMV8DB6uYZSIwywR8CySKSAfqfr+GldcY842VB2z+3oaxbytzIt/3GjnOrLZ+ZwGMMT8ZY1Zaw/nAeqBThdki9t1tLIWiE7Az6P0ujt3JwX5LoDKXiROR5SLyrYhcGIF8FYWb9xLrEPNdEelynMvWlrC3Z53O6w58FTS6rvdtdSr7PHW9X2ui4vfWAF+IyAoRuc6mTBWdLCLfi8inIjLAGldv962INCHwo/pe0Ghb96uIJANDgO8qTIrYd9d13CkbOBGZAqQBo4NGdzPG7BaRHsBXIrLWGLPVnoTlPgZmG2NKROR64BXgDJszVecK4F1jjC9oXH3ct1FHRMYSKBSnBo0+1dq3bYEvRWSD9Ze0XVYS+Pc+IiK/Aj4AetuYJxwXAEuMMcFHH7btVxFpRqBo3WSMyauLbULjOaLYDXQJet/ZGncUETkTuBOYYIwpKRtvjNlt/XcbsIBANY+kavMaY3KCMr4IDAt32Vp2PNu7ggqH8Dbs2+pU9nnqer+GTUQGEfgOTDTG5JSND9q3+4D3CZzisY0xJs8Yc8QangO4RaQ19XjfUvV3tk73q4i4CRSJ140x/wkxS+S+u3XZIGPXi8CR0zYCpz3KGssGVJhnCIEGtd4VxicBsdZwa2AzkW9oCydvh6Dhi4BvzS8NVz9auZOs4ZZ2ZrXm60egEVDs3LfWtpKpvMH1PI5uEFxmx349jrxdgS3AKRXGNwWaBw1/A5xjc9b2Zf/+BH5cd1j7OazvUF1mtaYnEGjHaFoP9qsAs4DHqpgnYt/dRnHqyRjjFZE/Ap8TuMJipjHmBxG5F1hujPkI+BfQDHhHRAB2GGMmAP2B50TET+AIbIYxZl09yPsnEZkAeAl8mTOsZQ+KyD+A/1qru9ccfdhsR1YI/GX2prG+uZY637ciMpvA1TetRWQX8HfAbX2WZ4E5BK4e2QIUAldb0+p0vx5H3ruAVsDT1vfWawK9h7YD3rfGuYA3jDGf2Zz1UuB3IuIFioArrO9DyO+QzVkh8AfYF8aYgqBF63y/WtKB3wBrRWS1Ne6vBP5QiPh3V7vwUEopVaXG0kahlFKqhrRQKKWUqpIWCqWUUlXSQqGUUqpKWiiUUkpVSQuFqvdE5EgY89xkdbdQW9u8UEROqsX1fXMCyx6x/ttRRN6tYr5EEfl9TbejVGW0UKiG4ibguAqFiDirmHwhUGuFwhhzSi2sY48x5tIqZkkEtFCoWqeFQkUN63kGC6xOEDeIyOtWH/x/AjoC80VkvjXveBFZKiIrReQdq4+csmcJPCQiK4HLRORaEfmv1VHdeyLSREROASYA/7KeOdBTRFKtjgvXiMj7IpJkrW+BBJ63sVxE1ovIcBH5j9Xv/31B2Y8EDd8mgecZfC8iM0J8zu5W9rUV1pEs1vMTRGSAiCyz8q0Rkd7ADKCnNe5fItJMAs9WWWmta2LQetaLyAsSeLbBFyISb03rJSJzrWwrRaSnNf5Waz+tEZF7avUfVtV/kb71XF/6OtEX1jMrCNxJe5hAXzUOYCmBDtog6BkBBLoD+Rqr6wXgNuCuoPn+ErTuVkHD9wE3WMOZwKVB09YAo63he7G6UiDQP9VD1vCNwB4Czw6IJdBLZ6sKn+FcAt0+lD375JiuFICPgKus4T8ELZuM1eUE8AQw2RqOAeKp0CUFgTuHWwTtky0EundIJnBHf6o17W1gijX8HXCRNRxH4ChtPPC8tawD+AQ43e7vhb7q7tUouvBQDcoyY8wuAKsrg2RgcYV5RhE4bbTE6mohhkBRKfNW0PBA66/2RAJduHxecYMikgAkGmMWWqNeAd4JmqWsm5K1wA/GmJ+s5bYR6IwtJ2jeM4GXjfXsExO6K4V04BJr+FXgoRDzLAXuFJHOwH+MMZutz3pUdOABETkd8BPoWrqdNe1HY8xqa3gFkCwizYFOxpj3rWzF1ucYT6BYrLLmb0ag11c7e6JVdUgLhYo2JUHDPkJ/h4XAg1qurGQdwX33ZAIXGmO+F5EMAkctNc3kr5DPX0m+cFTZt44x5g0R+Y5AR3BzJNDV/LYKs00G2gDDjDEeEckmcJQQnBkC+zG+is0J8KAx5rnjyK8aEG2jUA1FPoFHRELgSW/pItILQESaikifSpZrDvwkgS6cJ4danzHmMHBIRE6zpv0GWEjNfAlcXXaFlljPYa5gCYFOFKmQqZwEnt+xzRjzOPAhMIij9wEEej/dZxWJsUC3Y9f0CxN4ctousR4gJSKxVs7PgWlB7TydJPAsBtVIaKFQDcXzwGciMt8Ys59Ab7qzRWQNgdM0/SpZ7n8InJdfAmwIGv8mcKuIrLIadKcSaNxeA6QSaKc4bibQ0+hHwHLr1NktIWa7EfiDiKyl8ieRTQKyrHUMJPAIzBwCp9uyRORfwOtAmrWeqyp8vsr8hkDPxGsItKW0N8Z8AbwBLLXW9S5HFyTVwGnvsUoppaqkRxRKKaWqpIVCKaVUlbRQKKWUqpIWCqWUUlXSQqGUUqpKWiiUUkpVSQuFUkqpKv1/AXC3Tf2G0nMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "for value, bootstrap in conditions.items():\n",
    "    plt.plot(results[f'points_{bootstrap}'], results[f'energies_{bootstrap}'], label = f'{bootstrap}')\n",
    "plt.plot(results['points_np'], results['energies_np'], label = 'numpy')\n",
    "for condition in extrapolators.keys():\n",
    "    print(condition)\n",
    "    plt.plot(results[f'points_{condition}'], results[f'energies_{condition}'], label = condition)\n",
    "plt.legend()\n",
    "plt.title('Dissociation profile')\n",
    "plt.xlabel('Interatomic distance')\n",
    "plt.ylabel('Energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare number of evaluations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no bootstrapping\n",
      "Total evaluations for no bootstrapping:\n",
      "2931\n",
      "bootstrapping\n",
      "Total evaluations for bootstrapping:\n",
      "991\n",
      "np\n",
      "Total evaluations for np:\n",
      "0\n",
      "poly\n",
      "Total evaluations for poly:\n",
      "563\n",
      "diff_model\n",
      "Total evaluations for diff_model:\n",
      "1089\n"
     ]
    }
   ],
   "source": [
    "for condition, result_full in results_full.items():\n",
    "    print(condition)\n",
    "    print('Total evaluations for ' + condition + ':')\n",
    "    sum = 0\n",
    "    for key in results_full[condition].keys():\n",
    "        if condition is not 'np':\n",
    "                sum += results_full[condition][key]['raw_result']['cost_function_evals']\n",
    "        else:\n",
    "                sum = 0\n",
    "    print(sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>0.23.6</td></tr><tr><td>Terra</td><td>0.16.4</td></tr><tr><td>Aer</td><td>0.7.5</td></tr><tr><td>Ignis</td><td>0.5.2</td></tr><tr><td>Aqua</td><td>0.8.2</td></tr><tr><td>IBM Q Provider</td><td>0.11.1</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.8.6 | packaged by conda-forge | (default, Jan 25 2021, 23:21:18) \n",
       "[GCC 9.3.0]</td></tr><tr><td>OS</td><td>Linux</td></tr><tr><td>CPUs</td><td>8</td></tr><tr><td>Memory (Gb)</td><td>31.409000396728516</td></tr><tr><td colspan='2'>Tue Feb 23 18:07:35 2021 UTC</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2021.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
